EstuarineFishGillMicrobiota_part3_21.07.24
2024-07-21
We perform here different aspects of beta diversity analysis of the microbiota data gathered from fish gill mucus and bacterioplankton. Most relevant sections are the distance based 6.3.1 PCOA and 6.4.3 RDA analyses as well as the 6.3.5 Mantel test for global matrix correlations. Biomarker discovery via multipratt gives some indications for taxa specific for fish, bacterioplankton and seasonal sampling. Overall it seems adviced to run several biomarker discovery techniques in parallel but in this variable dataset they give little insides and are more meant as a backup for the network analysis.
1 Global Setup
1.1 Path
before you start: .rs.restartR()
1.2 Packages
1.3 Functions
1.4 Input Files
1.5 Tutorials
1.6 Setup Analysis
#-
6 Beta Diversity
6.1 PCA
CLR tranformed data in PCA analysis
##################
# SampleDist_PCAs# Publication
##################
#pslist <- pslistraw
for (i in names(pslist)[grepl("TSEc.l.r", names(pslist))]){
require(plyr)
require(ggrepel)
require(cowplot)
require(DESeq2)
require(matrixStats)
TSE <- pslist[[i]]
tryCatch({
g <- paste(i); print(i)
gg <- sub('TSEc.l.r_', '', g)
COL <- col.Palette$col.Palette.RepsSpecs[names(col.Palette$col.Palette.RepsSpecs) %in% TSE@colData$RepsSpecs]
A <- pcaplotRK_publication(TSE,intgroup = c("RepsSpecs"), pcX = 1, pcY = 2, title="", ellipse = F, ellipse.prob = 0.5, group_order = COL) +
scale_fill_manual(values = COL) +
scale_color_manual(values = COL) +
atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=6),
legend.text=element_text(size=6)) +
theme(
panel.grid.major = element_line(colour = "grey50"),
panel.grid.minor = element_line(colour = "grey50"))
SampleDist_PCA <- cowplot::plot_grid(A, labels = c(""), ncol = 1)
A <- plot_grid(SampleDist_PCA, ncol = 1)
plot(A)
ggsave(A, filename = paste(save_name, gg, "clrPCA.png"), path = pathPlots, device='png', dpi=300, width = 8,
height = 7)
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}## [1] "TSEc.l.r_SSU"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
## [1] "TSEc.l.r_WF"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
## [1] "TSEc.l.r_Fish"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
## [1] "TSEc.l.r_OE"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
## [1] "TSEc.l.r_GC"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
## [1] "TSEc.l.r_OEWF"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
## [1] "TSEc.l.r_GCWF"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."
6.2 CCA
###################
#CCA with relabund#
###################
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE|ps_Fish", names(pslist)[!grepl("WF", names(pslist))])]){
require(plyr)
require(ggrepel)
require(cowplot)
require(DESeq2)
require(matrixStats)
tryCatch({
ps <- pslist[[Dataset]]
Datasetname <- sub("ps_", "", Dataset)
if (Datasetname == "Fish") {
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% sample_data(ps)$Season]
GroupingVariable <- c("Season", "LOC")
} else {
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% sample_data(ps)$Season]
GroupingVariable <- c("Season", "LOC")
}
if (Datasetname %in% c("GC", "Fish")) {
ps <- prune_samples(sample_names(ps)[!grepl("GCSU21BBEB6", sample_names(ps))], ps)
} else {ps <- ps}
if (Datasetname == "GC") {
Traits <- c(
"Salinity","O2","SecchiDepth", "Temperature",
"Age", "FCF", "HSI", "SSI",
GroupingVariable)
print(Traits)
} else if (Datasetname %in% c("OE", "Fish")) {
Traits <- c(
"Salinity","O2","SecchiDepth", "Temperature",
"Cortison", GroupingVariable)
print(Traits)
} else {print("error")}
print(paste("Samples removed for low frequency below 10000 seqs in;", g, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
ps_Filter_relAbund <- ps_Filter %>%
transform_sample_counts(function(x) {x/sum(x)*100})
sample_data(ps_Filter)<-sample_data(ps_Filter)[ ,(names(sample_data(ps_Filter)) %in% Traits)]
'%!in%' <- function(x,y)!('%in%'(x,y))
print("NAs in TraitData")
print(table(is.na(sample_data(ps_Filter))))
include <- Traits[Traits != GroupingVariables]
'%!in%' <- function(x,y)!('%in%'(x,y))
Include <- Traits[Traits %!in% GroupingVariables]
df_Traits <- data.frame(sample_data(ps_Filter))
df_Traits <- df_Traits[names(df_Traits) %in% Traits]
df_Include <- df_Traits[names(df_Traits) %in% Include]
df_Include <- decostand(df_Include, method='standardize')
##################################################
#Determine which variables are useful for the plot#
##################################################
relAbund_table <- as.data.frame(otu_table(ps_Filter_relAbund, taxa_are_rows = FALSE))
##################################################
#Determine which variables are useful for the plot#
##################################################
abund_table.adonis <- vegan::adonis2(relAbund_table ~ ., data=df_Include)
print(Datasetname)
print(abund_table.adonis)
abund_table.adonis.sig <- abund_table.adonis[!is.na(abund_table.adonis$`Pr(>F)`),]
abund_table.adonis.sig <- abund_table.adonis.sig[abund_table.adonis.sig$`Pr(>F)` < 0.05,]
sig.variables <- row.names(abund_table.adonis.sig)
df_Significant <- df_Include[names(df_Include) %in% c(sig.variables)]
#########################
#Peform the CCA analysis#
#########################
perform_cca <- function(abund_table, meta_table) {
formula_text <- reformulate(response = "abund_table",
termlabels = names(meta_table))
tryCatch({
vegan.cca <- vegan::cca(formula_text, data = meta_table)
return(vegan.cca)
}, error = function(e) {
message("Error performing CCA: ", e$message)
return(NULL)
})
}
result.cca <- perform_cca(relAbund_table, df_Significant)
scores.cca<- vegan::scores(result.cca, display=c("sp","wa","lc","bp","cn"))
#########################
#Extract site data first#
df_sites <- data.frame(scores.cca$sites,df_Traits)
############
#Draw sites#
p <- ggplot(data=df_sites,aes(CCA1,CCA2,colour=Season, shape=LOC)) +
geom_point() +
scale_color_manual(values = COL) +
scale_fill_manual(values = COL) +
#ggrepel::geom_label_repel(mapping = aes_string(label = "LOC",
# fill = "LOC"), color = "white", show.legend = TRUE)
ggrepel::geom_text_repel(aes(label = LOC))
multiplier <- vegan:::ordiArrowMul(scores.cca$biplot)
df_arrows<- scores.cca$biplot*multiplier
df_arrows=as.data.frame(df_arrows)
p <- p + geom_segment(data=df_arrows,inherit.aes=F, mapping=(aes(x = 0, y = 0, xend = CCA1, yend = CCA2)), arrow = arrow(length = unit(0.2, "cm")),color="grey60",alpha=0.8) +
geom_text_repel(data=as.data.frame(df_arrows*1.2), inherit.aes=F, mapping=aes(CCA1, CCA2, label =
rownames(df_arrows)),color="grey30",alpha=0.99, size = 6)
p <- p + atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=6),
legend.text=element_text(size=6)) +
theme(
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "white"))
p <-p+ xlab(paste(sub("CAP1", "", ord_labels(result.cca)[1]))) +
ylab(paste(sub("CAP2", "", ord_labels(result.cca)[2])))
CCA_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
CCA_plot <- plot_grid(CCA_plot, ncol = 1)
plot(CCA_plot)
ggsave(CCA_plot, filename = paste(paste(save_name, "relabund_sigvariable_CCA", Datasetname, sep="_"), ".png", sep=""), path = pathPlots, device='png', dpi=300, width = 8,
height = 7)
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}## ERROR : error in evaluating the argument 'table' in selecting a method for function '%in%': could not find function "sample_data"
## ERROR : error in evaluating the argument 'table' in selecting a method for function '%in%': could not find function "sample_data"
## ERROR : error in evaluating the argument 'table' in selecting a method for function '%in%': could not find function "sample_data"6.3 Bray-Dist
We create a average distance dataset here by permutated subsampling
Bray_list <- list()
for (Dataset in c(names(pslist)[grepl("ps_GC|ps_OE|ps_Fish|ps_SSU", names(pslist))], "ps_WF")) {
require(phyloseq)
require(vegan)
require(plyr)
require(dplyr)
require(ggrepel)
require(cowplot)
require(ggordiplots)
Datasetname <- sub("ps_", "", Dataset)
#ps <- pslist[[Dataset]]
Samples <- sample_names(pslist[[Dataset]])
ps <- merge_phyloseq(pslist$ps_OE, pslist$ps_GC, pslist$ps_WF)
ps <- prune_samples(sample_names(ps) %in% Samples, ps)
Bray_list_length <- length(Bray_list)
if (Datasetname %in% c("GC", "OE", "GCWF", "OEWF", "Fish")) {
print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
} else if (Datasetname %in% c("SSU")) {
print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
} else if (Datasetname %in% c("WF")) {
print(paste("No Waterfiler Samples removed ;", Dataset, sep = ""))
print(rowSums(otu_table(ps)))
ps_Filter <- ps
}
ps_Filter_relAbund <- ps_Filter %>%
transform_sample_counts(function(x) {x/sum(x)*100})
###############
#RelAbund Dist#
###############
relAbund_dist_matrix <-
as.data.frame(otu_table(ps_Filter_relAbund, taxa_are_rows = FALSE)) %>%
vegdist(., method="bray")
##############
#Average Dist#
##############
smallest_group <- min(rowSums(otu_table(ps_Filter)))
print("Smallest group for avgdist")
print(smallest_group)
Avg_dist_matrix <- avgdist(as.data.frame(otu_table(ps_Filter,
taxa_are_rows = FALSE)), dmethod = "bray", sample = smallest_group)
relAbund_dist_dist_tibble <- relAbund_dist_matrix %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample) %>%
filter(name < sample )
Avg_dist_dist_tibble <- Avg_dist_matrix %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample) %>%
filter(name < sample)
Bray_list[[Bray_list_length+1]] <- relAbund_dist_matrix
names(Bray_list)[[Bray_list_length+1]] <- paste(Datasetname, "relAbund_dist_matrix", sep="_")
Bray_list[[Bray_list_length+2]] <- Avg_dist_matrix
names(Bray_list)[[Bray_list_length+2]] <- paste(Datasetname, "Avg_dist_matrix", sep="_")
Bray_list[[Bray_list_length+3]] <- ps_Filter
names(Bray_list)[[Bray_list_length+3]] <- paste("ps_Filter", Datasetname, sep="_")
}## Loading required package: phyloseq##
## Attaching package: 'phyloseq'## The following object is masked from 'package:SummarizedExperiment':
##
## distance## The following object is masked from 'package:Biobase':
##
## sampleNames## The following object is masked from 'package:GenomicRanges':
##
## distance## The following object is masked from 'package:IRanges':
##
## distance## Loading required package: vegan## Warning: package 'vegan' was built under R version 4.3.3## Loading required package: permute## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.3.1## This is vegan 2.6-6.1## Loading required package: ggordiplots## Warning: package 'ggordiplots' was built under R version 4.3.1## Loading required package: glue## Warning: package 'glue' was built under R version 4.3.1##
## Attaching package: 'glue'## The following object is masked from 'package:SummarizedExperiment':
##
## trim## The following object is masked from 'package:GenomicRanges':
##
## trim## The following object is masked from 'package:IRanges':
##
## trim## [1] "Samples removed for low frequency below 5000 seqs in;ps_SSU"
## GCAU21TWFL2 GCAU21MLEB1 WFSU21SSFL WFWI22MGEB
## 183 184 229 232
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "Samples removed for low frequency below 5000 seqs in;ps_Fish"
## GCAU21TWFL2 GCAU21MLEB1
## 183 184
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OE"
## named integer(0)
## [1] "Smallest group for avgdist"
## [1] 9071
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GC"
## GCAU21TWFL2 GCAU21MLEB1
## 45 46
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OEWF"
## WFSU21SSFL WFWI22MGEB
## 141 144
## [1] "Smallest group for avgdist"
## [1] 7739
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GCWF"
## GCAU21TWFL2 GCAU21MLEB1 WFSU21SSFL WFWI22MGEB
## 45 46 91 94
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "No Waterfiler Samples removed ;ps_WF"
## WFSU21MGEB WFSU21BBEB WFSU21SSFL WFSU21MLFL WFAU21TWEB WFWI22MGEB WFWI22MGFL
## 8983 8409 4572 17204 13785 3576 13431
## WFWI22BBEB WFWI22BBFL WFWI22MLFL WFWI22SSFL WFSP22MGEB WFSP22MGFL WFSP22BBEB
## 7739 26491 9265 24818 11908 31196 28110
## WFSP22BBFL WFSP22SSEB WFSP22SSFL WFSP22MLEB WFSP22MLFL
## 22786 17448 21514 9152 20503
## [1] "Smallest group for avgdist"
## [1] 3576pslist$BrayCurtis <- Bray_list6.3.1 PCOA
PCOA_list <- list()
for (Dataset in names(pslist$BrayCurtis)[grepl(
"OE_Avg_dist_matrix|GC_Avg_dist_matrix|SSU_Avg_dist_matrix|OE_relAbund_dist_matrix|GC_relAbund_dist_matrix|SSU_relAbund_dist_matrix|OEWF_Avg_dist_matrix|GCWF_Avg_dist_matrix", names(pslist$BrayCurtis))]) {
require(vegan)
require(plyr)
require(dplyr)
require(ggrepel)
require(cowplot)
Datasetname <- sub("_Avg_dist_matrix|_relAbund_dist_matrix", "", Dataset)
if (grepl("Avg_dist_matrix", Dataset)) {
BrayKind <- "avgdist"
} else if (grepl("relAbund_dist_matrix", Dataset)) {
BrayKind <- "relabund_dist"
}
Bray_list_species <- Bray_list[[Dataset]]
PCOA_list_length <- length(PCOA_list)
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
labels(Bray_list_species)]
set.seed(1)
pcoa <- cmdscale(Bray_list_species, eig=TRUE, add=TRUE)
positions <- pcoa$points
colnames(positions) <- c("pcoa1", "pcoa2")
percent_explained <- 100 * pcoa$eig / sum(pcoa$eig)
pretty_pe <- format(round(percent_explained, digits =1), nsmall=1, trim=TRUE)
labels <- c(glue("PCo Axis 1 ({pretty_pe[1]}%)"),
glue("PCo Axis 2 ({pretty_pe[2]}%)"))
if (Datasetname == "SSU") {
# p <- positions %>%
# as_tibble(rownames = "SampleID") %>%
# left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
# #left_join(data.frame(sample_data(Bray_list$ps_Filter_OE)), by = c("SampleID" = "SampleID" )) %>%
# ggplot(., aes(x=pcoa1, y=pcoa2, color=fSeason, shape=fSpecies)) +
# geom_point(size=4) +
# labs(x=labels[1], y=labels[2]) +
# #theme_bw() +
# theme(legend.position = "right") +
# guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
# shape = guide_legend(title.position="top", title.hjust = 0, title = "Location")) +
# scale_color_manual(values = col.Palette$col.Palette.Season) +
# scale_fill_manual(values = col.Palette$col.Palette.Season) +
# ggrepel::geom_text_repel(aes(label = fSpecies)) +
# coord_fixed(sqrt(percent_explained[2] / percent_explained[1]))
p <- positions %>%
as_tibble(rownames = "SampleID") %>%
left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
ggplot(aes(x=pcoa1, y=pcoa2, shape=fSpecies)) +
geom_point(aes(color=fSeason), size=4) +
labs(x=labels[1], y=labels[2]) +
theme(legend.position = "right") +
guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
shape = guide_legend(title.position="top", title.hjust = 0, title = "Species")) +
ggrepel::geom_text_repel(aes(label = fSpecies)) +
coord_fixed(sqrt(percent_explained[2] / percent_explained[1])) +
stat_ellipse(aes(fill=fSpecies), geom="polygon", alpha=0.25, colour = "black") +
scale_color_manual(values = col.Palette$col.Palette.Season) +
scale_fill_manual(values = col.Palette$col.Palette.Species)
} else if (Datasetname %in% c("GCWF", "OEWF")) {
p <- positions %>%
as_tibble(rownames = "SampleID") %>%
left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
ggplot(., aes(x=pcoa1, y=pcoa2)) +
geom_point(size=4, aes(color=fSeason, shape=fLOC)) +
labs(x=labels[1], y=labels[2]) +
theme(legend.position = "right") +
#guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
# shape = guide_legend(title.position="top", title.hjust = 0, title = "Species")) +
#ggrepel::geom_text_repel(aes(label = SampleID)) +
coord_fixed(sqrt(percent_explained[2] / percent_explained[1])) +
stat_ellipse(aes(fill=fSpecies), geom="polygon", alpha=0.25, colour = "black") +
scale_color_manual(values = col.Palette$col.Palette.Season) +
scale_fill_manual(values = col.Palette$col.Palette.Species) +
scale_shape_manual(values=c(15, 16, 17, 8, 25)) +
#ggrepel::geom_text_repel(aes(label = fSpecies))+
coord_fixed(sqrt(percent_explained[2] / percent_explained[1]))
} else {
p <- positions %>%
as_tibble(rownames = "SampleID") %>%
left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
#left_join(data.frame(sample_data(Bray_list$ps_Filter_OE)), by = c("SampleID" = "SampleID" )) %>%
ggplot(., aes(x=pcoa1, y=pcoa2, color=fSeason, shape=fLOC)) +
geom_point(size=4) +
labs(x=labels[1], y=labels[2]) +
theme_bw()+ theme(legend.position = "right") +
guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
shape = guide_legend(title.position="top", title.hjust = 0, title = "Location")) +
scale_color_manual(values = col.Palette$col.Palette.Season) +
scale_fill_manual(values = col.Palette$col.Palette.Season) +
scale_shape_manual(values=c(15, 16, 17, 8, 25)) +
ggrepel::geom_text_repel(aes(label = fLOC))+
coord_fixed(sqrt(percent_explained[2] / percent_explained[1]))
}
#p <- p + annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1, size = 5,
# label = paste("Stress =", round(stress, 3)),
# color = c("darkred", "darkred", "darkred"), fontface = "bold")
#p + theme(panel.grid = element_blank())
p <- p + atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=6),
legend.text=element_text(size=6)) +
theme(
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "white"))
p <- p + theme(legend.position = "bottom", legend.box = "horizontal")
PCOA_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
PCOA_plot <- plot_grid(PCOA_plot, ncol = 1)
#plot(PCOA_plot)
ggsave(PCOA_plot, filename = paste(paste(save_name, BrayKind, "PCOA", Datasetname, sep="_"), ".png", sep=""),
path = pathPlots, device='png', dpi=300, width = 8,height = 7)
PCOA_list[[PCOA_list_length+1]] <- PCOA_plot
names(PCOA_list)[[PCOA_list_length+1]] <- paste("PCOA",BrayKind, Datasetname, sep="_")
}## Warning: ggrepel: 191 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 195 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 69 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 70 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 34 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 32 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.PCOA_list$PCOA_avgdist_SSU## Warning: ggrepel: 212 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
6.3.2 NMDS
NMDS_list <- list()
for (Dataset in names(pslist$BrayCurtis)[grepl(
"OE_Avg_dist_matrix|GC_Avg_dist_matrix|SSU_Avg_dist_matrix|OE_relAbund_dist_matrix|GC_relAbund_dist_matrix|SSU_relAbund_dist_matrix", names(pslist$BrayCurtis))]) {
require(vegan)
require(plyr)
require(dplyr)
require(ggrepel)
require(cowplot)
Datasetname <- sub("_Avg_dist_matrix|_relAbund_dist_matrix", "", Dataset)
if (grepl("Avg_dist_matrix", Dataset)) {
BrayKind <- "avgdist"
} else if (grepl("relAbund_dist_matrix", Dataset)) {
BrayKind <- "relabund_dist"
}
Bray_list_species <- Bray_list[[Dataset]]
NMDS_list_length <- length(NMDS_list)
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
labels(Bray_list_species)]
set.seed(1)
nmds <- metaMDS(Bray_list_species)
stress <- nmds$stress
if (Datasetname == "SSU") {
p <- scores(nmds) %>%
as_tibble(rownames = "samples") %>%
left_join(SAMDF16S, by = c("samples" = "SampleID" )) %>%
ggplot(aes(x=NMDS1, y=NMDS2, color=fSeason, shape=fLOC)) +
geom_point(size=4)+
scale_color_manual(values = col.Palette$col.Palette.Season) +
scale_fill_manual(values = col.Palette$col.Palette.Season) +
ggrepel::geom_text_repel(aes(label = fSpecies))
} else {
p <- scores(nmds) %>%
as_tibble(rownames = "samples") %>%
left_join(SAMDF16S, by = c("samples" = "SampleID" )) %>%
ggplot(aes(x=NMDS1, y=NMDS2, color=fSeason, shape=fLOC)) +
geom_point(size=4)+
scale_color_manual(values = col.Palette$col.Palette.Season) +
scale_fill_manual(values = col.Palette$col.Palette.Season) +
scale_shape_manual(values=c(15, 16, 17, 8, 25)) +
ggrepel::geom_text_repel(aes(label = fLOC))
}
p <- p + annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1, size = 5,
label = paste("Stress =", round(stress, 3)),
color = c("darkred", "darkred", "darkred"), fontface = "bold")
p <- p + atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=6),
legend.text=element_text(size=6)) +
theme(
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "white"))
NMDS_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
NMDS_plot <- plot_grid(NMDS_plot, ncol = 1)
#plot(NMDS_plot)
ggsave(NMDS_plot, filename = paste(paste(save_name, BrayKind, "NMDS", Datasetname, sep="_"), ".png", sep=""),
path = pathPlots, device='png', dpi=300, width = 8,height = 7)
NMDS_list[[NMDS_list_length+1]] <- NMDS_plot
names(NMDS_list)[[NMDS_list_length+1]] <- paste("NMDS",BrayKind, Datasetname, sep="_")
}## Run 0 stress 0.1598705
## Run 1 stress 0.1589304
## ... New best solution
## ... Procrustes: rmse 0.01986886 max resid 0.2209935
## Run 2 stress 0.1603834
## Run 3 stress 0.1673076
## Run 4 stress 0.157282
## ... New best solution
## ... Procrustes: rmse 0.03117025 max resid 0.2195268
## Run 5 stress 0.1872517
## Run 6 stress 0.1598648
## Run 7 stress 0.1587316
## Run 8 stress 0.1568198
## ... New best solution
## ... Procrustes: rmse 0.01454307 max resid 0.178522
## Run 9 stress 0.1632134
## Run 10 stress 0.1577244
## Run 11 stress 0.1629727
## Run 12 stress 0.162905
## Run 13 stress 0.1625196
## Run 14 stress 0.1646227
## Run 15 stress 0.1600792
## Run 16 stress 0.1625054
## Run 17 stress 0.1604759
## Run 18 stress 0.1706289
## Run 19 stress 0.1629655
## Run 20 stress 0.1606143
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 3: no. of iterations >= maxit
## 17: stress ratio > sratmax## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 156 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Run 0 stress 0.158884
## Run 1 stress 0.157863
## ... New best solution
## ... Procrustes: rmse 0.0201235 max resid 0.2197741
## Run 2 stress 0.1598936
## Run 3 stress 0.1662939
## Run 4 stress 0.1573704
## ... New best solution
## ... Procrustes: rmse 0.02908927 max resid 0.2207995
## Run 5 stress 0.4174065
## Run 6 stress 0.1560476
## ... New best solution
## ... Procrustes: rmse 0.02829525 max resid 0.2234037
## Run 7 stress 0.1566534
## Run 8 stress 0.1579184
## Run 9 stress 0.1622481
## Run 10 stress 0.1569373
## Run 11 stress 0.1642362
## Run 12 stress 0.1629476
## Run 13 stress 0.1613905
## Run 14 stress 0.1640273
## Run 15 stress 0.159106
## Run 16 stress 0.1614624
## Run 17 stress 0.1594659
## Run 18 stress 0.1681732
## Run 19 stress 0.1615566
## Run 20 stress 0.1595685
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 2: no. of iterations >= maxit
## 15: stress ratio > sratmax
## 3: scale factor of the gradient < sfgrmin## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 162 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Run 0 stress 0.1406636
## Run 1 stress 0.1642062
## Run 2 stress 0.1614026
## Run 3 stress 0.1697188
## Run 4 stress 0.1544251
## Run 5 stress 0.1627928
## Run 6 stress 0.1698772
## Run 7 stress 0.1398245
## ... New best solution
## ... Procrustes: rmse 0.01506232 max resid 0.1027648
## Run 8 stress 0.1541418
## Run 9 stress 0.1609619
## Run 10 stress 0.1570996
## Run 11 stress 0.1423337
## Run 12 stress 0.1454004
## Run 13 stress 0.1562743
## Run 14 stress 0.157439
## Run 15 stress 0.1666092
## Run 16 stress 0.1434355
## Run 17 stress 0.1397957
## ... New best solution
## ... Procrustes: rmse 0.001201193 max resid 0.01112539
## Run 18 stress 0.1738992
## Run 19 stress 0.1426277
## Run 20 stress 0.1657134
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 13: stress ratio > sratmax
## 7: scale factor of the gradient < sfgrmin## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 74 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Run 0 stress 0.1400518
## Run 1 stress 0.1623834
## Run 2 stress 0.1605131
## Run 3 stress 0.1687144
## Run 4 stress 0.1530577
## Run 5 stress 0.1626125
## Run 6 stress 0.1692477
## Run 7 stress 0.1391221
## ... New best solution
## ... Procrustes: rmse 0.01506125 max resid 0.1005011
## Run 8 stress 0.1536073
## Run 9 stress 0.1598388
## Run 10 stress 0.1560175
## Run 11 stress 0.1414518
## Run 12 stress 0.14463
## Run 13 stress 0.1556535
## Run 14 stress 0.1575412
## Run 15 stress 0.1662978
## Run 16 stress 0.1426365
## Run 17 stress 0.1391944
## ... Procrustes: rmse 0.001685777 max resid 0.01933529
## Run 18 stress 0.1720325
## Run 19 stress 0.141849
## Run 20 stress 0.1606002
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 16: stress ratio > sratmax
## 4: scale factor of the gradient < sfgrmin## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 73 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Run 0 stress 0.1636024
## Run 1 stress 0.1546993
## ... New best solution
## ... Procrustes: rmse 0.07776629 max resid 0.6089917
## Run 2 stress 0.1672091
## Run 3 stress 0.1602318
## Run 4 stress 0.1589382
## Run 5 stress 0.156329
## Run 6 stress 0.1682938
## Run 7 stress 0.1715344
## Run 8 stress 0.1740996
## Run 9 stress 0.171455
## Run 10 stress 0.1590428
## Run 11 stress 0.1673151
## Run 12 stress 0.1619349
## Run 13 stress 0.1730834
## Run 14 stress 0.1614785
## Run 15 stress 0.1594615
## Run 16 stress 0.154597
## ... New best solution
## ... Procrustes: rmse 0.01151229 max resid 0.09751585
## Run 17 stress 0.155883
## Run 18 stress 0.1698926
## Run 19 stress 0.1602295
## Run 20 stress 0.1851819
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 5: no. of iterations >= maxit
## 15: stress ratio > sratmax## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: ggrepel: 51 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Run 0 stress 0.1638112
## Run 1 stress 0.1549974
## ... New best solution
## ... Procrustes: rmse 0.07744348 max resid 0.6090912
## Run 2 stress 0.1686948
## Run 3 stress 0.1623907
## Run 4 stress 0.1586151
## Run 5 stress 0.1630234
## Run 6 stress 0.172169
## Run 7 stress 0.1721439
## Run 8 stress 0.1677825
## Run 9 stress 0.1682698
## Run 10 stress 0.1653278
## Run 11 stress 0.1675991
## Run 12 stress 0.161257
## Run 13 stress 0.1747388
## Run 14 stress 0.1596532
## Run 15 stress 0.1672545
## Run 16 stress 0.1683737
## Run 17 stress 0.1545751
## ... New best solution
## ... Procrustes: rmse 0.01840541 max resid 0.1318633
## Run 18 stress 0.1703611
## Run 19 stress 0.1622102
## Run 20 stress 0.1855765
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 16: stress ratio > sratmax## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## ggrepel: 51 unlabeled data points (too many overlaps). Consider increasing max.overlaps6.3.3 PERMANOVA
https://www.nicholas-ollberding.com/post/introduction-to-the-statistical-analysis-of-microbiome-data-in-r/ While PCA is an exploratory data visualization tool, we can test whether the samples cluster beyond that expected by sampling variability using permutational multivariate analysis of variance (PERMANOVA). It does this by partitioning the sums of squares for the within- and between-cluster components using the concept of centroids. Many permutations of the data (i.e. random shuffling) are used to generate the null distribution. The test from ADONIS can be confounded by differences in dispersion (or spread)…so we want to check this as well.
#########################################################
#Test null hypotheses for non-rarified and rarified data#
#########################################################
# treatment_sample <- sample(treatment_size)
# adonis2(Bray_list[[paste(Datasetname, "relAbund_dist_matrix", sep="_")]]~ treatment_sample)$`Pr(>F)`[1]
# #0.269
# adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]~ treatment_sample)$`Pr(>F)`[1]
# #0.275
#####################
#Fish vs Waterfilter#
#####################
PERMANOVA_list <- list()
for (Dataset in names(pslist)[grepl("ps_OEWF|ps_GCWF", names(pslist))]){
require(vegan)
PERMANOVA_list[[Dataset]] <- list()
Datasetname <- sub("ps_", "", Dataset)
Bray_list <- pslist$BrayCurtis
#Other possibility#
#Generate distance matrix
#ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
#Generate distance matrix
#clr_dist_matrix <- phyloseq::distance(ps_clr, method = "euclidean")
#ADONIS test
print(Datasetname)
ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Kind)
print(ADONIS)
#Dispersion test and plot
dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Kind)
dispr
plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
boxplot(dispr, main = "", xlab = "")
vegan::permutest(dispr)
BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Kind))
#Cant trust Adonis for both
if (BetaDisper$`Pr(>F)`[1] > 0.05) {
print(Datasetname)
print("Trust Adonis")
} else {
print(Datasetname)
print("Can not trust Adonis")
}
print(BetaDisper )
print("Next analysis")
PERMANOVA_list[[Dataset]][[paste("Kind")]]<- ADONIS
}## [1] "OEWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 1
## Residual 153
## Total 154
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 6.813
## Residual 36.228
## Total 43.041
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 0.15829
## Residual 0.84171
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 28.773
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "OEWF"
## [1] "Trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.0022 0.0021677 0.0768 0.7821
## Residuals 153 4.3186 0.0282261
## [1] "Next analysis"
## [1] "GCWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 1
## Residual 101
## Total 102
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 7.0305
## Residual 20.5523
## Total 27.5828
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 0.25489
## Residual 0.74511
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 34.55
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "GCWF"
## [1] "Trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.01827 0.018273 0.7806 0.379
## Residuals 101 2.36418 0.023408
## [1] "Next analysis" ##############################
#Fish Seasonal vs Waterfilter#
##############################
for (Dataset in names(pslist)[grepl("ps_OEWF|ps_GCWF", names(pslist))]){
require(vegan)
PERMANOVA_list[[paste(Dataset, Season, sep="_")]] <- list()
Datasetname <- sub("ps_", "", Dataset)
#Other possibility#
#Generate distance matrix
#ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
#Generate distance matrix
#clr_dist_matrix <- phyloseq::distance(ps_clr, method = "euclidean")
#ADONIS test
metadata <- sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])
print(Datasetname)
ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ gsub(" \\d+$", "", metadata$Reps))
print(ADONIS)
#Dispersion test and plot
dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], gsub(" \\d+$", "", metadata$Reps))
dispr
plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
boxplot(dispr, main = "", xlab = "")
vegan::permutest(dispr)
BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], gsub(" \\d+$", "", metadata$Reps)))
#Cant trust Adonis for both
if (BetaDisper$`Pr(>F)`[1] > 0.05) {
print(Datasetname)
print("Trust Adonis")
} else {
print(Datasetname)
print("Can not trust Adonis")
}
print(BetaDisper )
print("Next analysis")
PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("ADONIS")]] <- list()
PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("ADONIS")]]<- ADONIS
print(Datasetname)
metadata$RepsSeason <- gsub(" \\d+$", "", metadata$Reps)
dist_matrix <- Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]
PairwiseADONIS <- pairwiseAdonis::pairwise.adonis2(dist_matrix~RepsSeason, data=data.frame(metadata))
PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("Pairwise")]] <- list()
PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("Pairwise")]]<- PairwiseADONIS
PairwiseADONIS_df <- list()
for (i in names(PairwiseADONIS[-1])) {
PairwiseADONIS_df[[i]][1] <- PairwiseADONIS[[i]]$`Pr(>F)`[1]
PairwiseADONIS_df[[i]][2] <- PairwiseADONIS[[i]]$`F`[1]
#names(PairwiseADONIS_df)[[i]] <- i
}
PairwiseADONIS_df <- list()
for (i in names(PairwiseADONIS)[-1]) {
PairwiseADONIS_df[[i]] <- data.frame(
p_value = PairwiseADONIS[[i]]$`Pr(>F)`[1],
F_statistic = PairwiseADONIS[[i]]$`F`[1]
)
}
PairwiseADONIS_df <- do.call(rbind, PairwiseADONIS_df)
print(PairwiseADONIS_df)
}## [1] "OEWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ gsub(" \\d+$", "", metadata$Reps))
## Df SumOfSqs R2 F Pr(>F)
## gsub(" \\\\d+$", "", metadata$Reps) 4 18.714 0.43478 28.846 0.001 ***
## Residual 150 24.328 0.56522
## Total 154 43.041 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "OEWF"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 4 1.0837 0.270915 24.721 9.904e-16 ***
## Residuals 150 1.6439 0.010959
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "OEWF"
## p_value F_statistic
## SU_vs_AU 0.001 26.214822
## SU_vs_WI 0.001 14.329574
## SU_vs_SP 0.001 6.215647
## SU_vs_WF 0.001 34.369533
## AU_vs_WI 0.001 34.717334
## AU_vs_SP 0.001 39.381833
## AU_vs_WF 0.001 18.085671
## WI_vs_SP 0.001 19.891102
## WI_vs_WF 0.001 47.848857
## SP_vs_WF 0.001 51.100772
## [1] "GCWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ gsub(" \\d+$", "", metadata$Reps))
## Df SumOfSqs R2 F Pr(>F)
## gsub(" \\\\d+$", "", metadata$Reps) 4 9.7506 0.3535 13.396 0.001 ***
## Residual 98 17.8322 0.6465
## Total 102 27.5828 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "GCWF"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 4 0.46918 0.117296 6.0819 0.0002054 ***
## Residuals 98 1.89002 0.019286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "GCWF"
## p_value F_statistic
## SU_vs_AU 0.001 6.165993
## SU_vs_WI 0.001 3.731333
## SU_vs_SP 0.001 7.507221
## SU_vs_WF 0.001 20.103508
## AU_vs_WI 0.001 4.375215
## AU_vs_SP 0.001 3.805386
## AU_vs_WF 0.001 29.802728
## WI_vs_SP 0.001 4.289426
## WI_vs_WF 0.001 16.873265
## SP_vs_WF 0.001 36.695406 df <- PERMANOVA_list$ps_OEWF_SU_AU_WI_SP$Pairwise
PairwiseADONIS_df <- list()
for (i in names(df)[-1]) {
PairwiseADONIS_df[[i]] <- data.frame(
p_value = df[[i]]$`Pr(>F)`[1],
F_statistic = df[[i]]$`F`[1]
)
}
PairwiseADONIS_df <- do.call(rbind, PairwiseADONIS_df)
PairwiseADONIS_df[grepl("WF", rownames(PairwiseADONIS_df)),]## p_value F_statistic
## SU_vs_WF 0.001 34.36953
## AU_vs_WF 0.001 18.08567
## WI_vs_WF 0.001 47.84886
## SP_vs_WF 0.001 51.10077 df <- PERMANOVA_list$ps_GCWF_SU_AU_WI_SP$Pairwise
PairwiseADONIS_df <- list()
for (i in names(df)[-1]) {
PairwiseADONIS_df[[i]] <- data.frame(
p_value = df[[i]]$`Pr(>F)`[1],
F_statistic = df[[i]]$`F`[1]
)
}
PairwiseADONIS_df <- do.call(rbind, PairwiseADONIS_df)
PairwiseADONIS_df[grepl("WF", rownames(PairwiseADONIS_df)),]## p_value F_statistic
## SU_vs_WF 0.001 20.10351
## AU_vs_WF 0.001 29.80273
## WI_vs_WF 0.001 16.87326
## SP_vs_WF 0.001 36.69541 ##################
#For Fish vs Fish#
##################
#Pairwise Test
for (Dataset in names(pslist)[grepl("ps_Fish", names(pslist))]){
require(vegan)
Datasetname <- sub("ps_", "", Dataset)
#Other Approach
#Generate distance matrix
# ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
# set.seed(123)
# clr_dist_matrix <- phyloseq::distance(ps_clr, method ="euclidean")
print(Datasetname)
ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Species)
print(ADONIS)
dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Species)
dispr
plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
boxplot(dispr, main = "", xlab = "")
vegan::permutest(dispr)
BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Species))
#Cant trust Adonis for both
if (BetaDisper$`Pr(>F)`[1] > 0.05) {
print(Datasetname)
print("Trust Adonis")
} else {
print(Datasetname)
print("Can not trust Adonis")
}
print(BetaDisper )
print("Next analysis")
PERMANOVA_list[[Dataset]][[paste("Species")]]<- ADONIS
}## [1] "Fish"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 1
## Residual 222
## Total 223
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 1.353
## Residual 49.747
## Total 51.100
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 0.02647
## Residual 0.97353
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 6.0362
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Fish"
## [1] "Trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.1038 0.10381 3.6425 0.05761 .
## Residuals 222 6.3270 0.02850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis" ###################
#For Fish Seasonal#
###################
#Pairwise Test
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_OE|ps_GC", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")){
require(vegan)
PERMANOVA_list[[Dataset]][[paste("Season")]] <- list()
Datasetname <- sub("ps_", "", Dataset)
#Other Approach
#Generate distance matrix
# ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
# set.seed(123)
# clr_dist_matrix <- phyloseq::distance(ps_clr, method ="euclidean")
print(Datasetname)
ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Season)
print(ADONIS)
dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Season)
dispr
plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
boxplot(dispr, main = "", xlab = "")
vegan::permutest(dispr)
BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Season))
#Cant trust Adonis for both
if (BetaDisper$`Pr(>F)`[1] > 0.05) {
print(Datasetname)
print("Trust Adonis")
} else {
print(Datasetname)
print("Can not trust Adonis")
}
print(BetaDisper )
print("Next analysis")
PERMANOVA_list[[Dataset]][[paste("Season")]][[paste("Adonis")]] <- list()
PERMANOVA_list[[Dataset]][[paste("Season")]][[paste("Adonis")]] <- ADONIS
metadata <- sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])
dist_matrix <- Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]
print(Datasetname)
PairwiseADONIS <- pairwiseAdonis::pairwise.adonis2(dist_matrix~Season,data=data.frame(metadata))
PairwiseADONIS
PERMANOVA_list[[Dataset]][[paste("Season")]][[paste("Pairwise")]]<-PairwiseADONIS
}## [1] "OE"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 3
## Residual 134
## Total 137
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 11.909
## Residual 20.664
## Total 32.573
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.36562
## Residual 0.63438
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 25.743
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "OE"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.98489 0.32830 29.441 1.097e-14 ***
## Residuals 134 1.49422 0.01115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "OE"
## [1] "GC"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 3
## Residual 82
## Total 85
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 2.7221
## Residual 14.2016
## Total 16.9237
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.16084
## Residual 0.83916
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 5.2391
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "GC"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.38924 0.129746 6.1046 0.0008426 ***
## Residuals 82 1.74281 0.021254
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "GC"
## [1] "WF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 3
## Residual 15
## Total 18
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 1.4831
## Residual 2.8284
## Total 4.3115
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.34398
## Residual 0.65602
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 2.6218
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "WF"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.16145 0.053816 4.2135 0.02386 *
## Residuals 15 0.19159 0.012772
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "WF" ##################
#For Fish Spatial#
##################
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_OE|ps_GC", names(pslist)[!grepl("WF", names(pslist))])]){
require(vegan)
PERMANOVA_list[[Dataset]][[paste("Replicates")]] <- list()
Datasetname <- sub("ps_", "", Dataset)
#Other Approach
#Generate distance matrix
# ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
# set.seed(123)
# clr_dist_matrix <- phyloseq::distance(ps_clr, method ="euclidean")
print(Datasetname)
ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~
sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$fReplicates)
print(ADONIS)
dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]],
sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$fReplicates)
dispr
plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
boxplot(dispr, main = "", xlab = "")
vegan::permutest(dispr)
BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]],
sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$fReplicates))
#Cant trust Adonis for both
if (BetaDisper$`Pr(>F)`[1] > 0.05) {
print(Datasetname)
print("Trust Adonis")
} else {
print(Datasetname)
print("Can not trust Adonis")
}
print(BetaDisper )
print("Next analysis")
PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Adonis")]] <- list()
PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Adonis")]] <- ADONIS
metadata <- sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])
dist_matrix <- Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]
PairwiseADONIS <- pairwiseAdonis::pairwise.adonis2(dist_matrix~fReplicates,data=data.frame(metadata))
PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Pairwise")]] <- list()
PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Pairwise")]]<- PairwiseADONIS
PairwiseADONIS_df <- list()
for (i in names(PairwiseADONIS[-1])) {
PairwiseADONIS_df[[i]] <- PairwiseADONIS[[i]]$`Pr(>F)`[1]
#names(PairwiseADONIS_df)[[i]] <- i
}
PairwiseADONIS_df <- PairwiseADONIS_df %>%
do.call(rbind, .) %>%
as.data.frame() %>%
tibble::rownames_to_column(var = "Comparison") %>%
dplyr::rename(p=2) %>%
filter(p < 0.05)
print(PairwiseADONIS_df)
}## [1] "OE"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 19
## Residual 118
## Total 137
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 19.874
## Residual 12.699
## Total 32.573
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 0.61015
## Residual 0.38985
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 9.7198
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "OE"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 19 0.76993 0.040523 4.652 7.784e-08 ***
## Residuals 118 1.02789 0.008711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## Comparison p
## 1 OESU21MG_vs_OESU21SS 0.001
## 2 OESU21MG_vs_OESU21TW 0.002
## 3 OESU21MG_vs_OESU21ML 0.003
## 4 OESU21MG_vs_OEAU21MG 0.001
## 5 OESU21MG_vs_OEAU21BB 0.001
## 6 OESU21MG_vs_OEAU21SS 0.003
## 7 OESU21MG_vs_OEAU21TW 0.001
## 8 OESU21MG_vs_OEAU21ML 0.001
## 9 OESU21MG_vs_OEWI22BB 0.003
## 10 OESU21MG_vs_OEWI22MG 0.001
## 11 OESU21MG_vs_OEWI22SS 0.001
## 12 OESU21MG_vs_OEWI22TW 0.002
## 13 OESU21MG_vs_OEWI22ML 0.003
## 14 OESU21MG_vs_OESP22MG 0.002
## 15 OESU21MG_vs_OESP22BB 0.002
## 16 OESU21MG_vs_OESP22SS 0.001
## 17 OESU21MG_vs_OESP22TW 0.003
## 18 OESU21MG_vs_OESP22ML 0.001
## 19 OESU21BB_vs_OESU21SS 0.017
## 20 OESU21BB_vs_OESU21TW 0.002
## 21 OESU21BB_vs_OESU21ML 0.001
## 22 OESU21BB_vs_OEAU21MG 0.001
## 23 OESU21BB_vs_OEAU21BB 0.001
## 24 OESU21BB_vs_OEAU21SS 0.001
## 25 OESU21BB_vs_OEAU21TW 0.001
## 26 OESU21BB_vs_OEAU21ML 0.003
## 27 OESU21BB_vs_OEWI22BB 0.004
## 28 OESU21BB_vs_OEWI22MG 0.001
## 29 OESU21BB_vs_OEWI22SS 0.001
## 30 OESU21BB_vs_OEWI22TW 0.002
## 31 OESU21BB_vs_OEWI22ML 0.002
## 32 OESU21BB_vs_OESP22MG 0.002
## 33 OESU21BB_vs_OESP22BB 0.038
## 34 OESU21BB_vs_OESP22SS 0.006
## 35 OESU21BB_vs_OESP22TW 0.001
## 36 OESU21BB_vs_OESP22ML 0.011
## 37 OESU21SS_vs_OESU21TW 0.001
## 38 OESU21SS_vs_OESU21ML 0.004
## 39 OESU21SS_vs_OEAU21MG 0.001
## 40 OESU21SS_vs_OEAU21BB 0.001
## 41 OESU21SS_vs_OEAU21SS 0.001
## 42 OESU21SS_vs_OEAU21TW 0.002
## 43 OESU21SS_vs_OEAU21ML 0.002
## 44 OESU21SS_vs_OEWI22BB 0.002
## 45 OESU21SS_vs_OEWI22MG 0.001
## 46 OESU21SS_vs_OEWI22SS 0.002
## 47 OESU21SS_vs_OEWI22TW 0.001
## 48 OESU21SS_vs_OEWI22ML 0.001
## 49 OESU21SS_vs_OESP22MG 0.001
## 50 OESU21SS_vs_OESP22BB 0.001
## 51 OESU21SS_vs_OESP22SS 0.001
## 52 OESU21SS_vs_OESP22TW 0.001
## 53 OESU21SS_vs_OESP22ML 0.008
## 54 OESU21TW_vs_OESU21ML 0.029
## 55 OESU21TW_vs_OEAU21MG 0.001
## 56 OESU21TW_vs_OEAU21BB 0.002
## 57 OESU21TW_vs_OEAU21SS 0.001
## 58 OESU21TW_vs_OEAU21TW 0.002
## 59 OESU21TW_vs_OEAU21ML 0.002
## 60 OESU21TW_vs_OEWI22BB 0.001
## 61 OESU21TW_vs_OEWI22MG 0.001
## 62 OESU21TW_vs_OEWI22SS 0.001
## 63 OESU21TW_vs_OEWI22TW 0.001
## 64 OESU21TW_vs_OEWI22ML 0.001
## 65 OESU21TW_vs_OESP22MG 0.002
## 66 OESU21TW_vs_OESP22BB 0.001
## 67 OESU21TW_vs_OESP22SS 0.001
## 68 OESU21TW_vs_OESP22TW 0.003
## 69 OESU21TW_vs_OESP22ML 0.001
## 70 OESU21ML_vs_OEAU21MG 0.002
## 71 OESU21ML_vs_OEAU21BB 0.002
## 72 OESU21ML_vs_OEAU21SS 0.003
## 73 OESU21ML_vs_OEAU21TW 0.003
## 74 OESU21ML_vs_OEAU21ML 0.002
## 75 OESU21ML_vs_OEWI22BB 0.002
## 76 OESU21ML_vs_OEWI22MG 0.001
## 77 OESU21ML_vs_OEWI22SS 0.001
## 78 OESU21ML_vs_OEWI22TW 0.004
## 79 OESU21ML_vs_OEWI22ML 0.002
## 80 OESU21ML_vs_OESP22MG 0.003
## 81 OESU21ML_vs_OESP22BB 0.002
## 82 OESU21ML_vs_OESP22SS 0.004
## 83 OESU21ML_vs_OESP22TW 0.002
## 84 OESU21ML_vs_OESP22ML 0.004
## 85 OEAU21MG_vs_OEAU21BB 0.001
## 86 OEAU21MG_vs_OEAU21SS 0.001
## 87 OEAU21MG_vs_OEAU21TW 0.002
## 88 OEAU21MG_vs_OEAU21ML 0.002
## 89 OEAU21MG_vs_OEWI22BB 0.001
## 90 OEAU21MG_vs_OEWI22MG 0.002
## 91 OEAU21MG_vs_OEWI22SS 0.001
## 92 OEAU21MG_vs_OEWI22TW 0.004
## 93 OEAU21MG_vs_OEWI22ML 0.002
## 94 OEAU21MG_vs_OESP22MG 0.001
## 95 OEAU21MG_vs_OESP22BB 0.001
## 96 OEAU21MG_vs_OESP22SS 0.006
## 97 OEAU21MG_vs_OESP22TW 0.001
## 98 OEAU21MG_vs_OESP22ML 0.004
## 99 OEAU21BB_vs_OEAU21SS 0.004
## 100 OEAU21BB_vs_OEAU21TW 0.001
## 101 OEAU21BB_vs_OEAU21ML 0.001
## 102 OEAU21BB_vs_OEWI22BB 0.001
## 103 OEAU21BB_vs_OEWI22MG 0.001
## 104 OEAU21BB_vs_OEWI22SS 0.001
## 105 OEAU21BB_vs_OEWI22TW 0.001
## 106 OEAU21BB_vs_OEWI22ML 0.001
## 107 OEAU21BB_vs_OESP22MG 0.002
## 108 OEAU21BB_vs_OESP22BB 0.003
## 109 OEAU21BB_vs_OESP22SS 0.001
## 110 OEAU21BB_vs_OESP22TW 0.002
## 111 OEAU21BB_vs_OESP22ML 0.001
## 112 OEAU21SS_vs_OEAU21TW 0.007
## 113 OEAU21SS_vs_OEAU21ML 0.001
## 114 OEAU21SS_vs_OEWI22BB 0.001
## 115 OEAU21SS_vs_OEWI22MG 0.003
## 116 OEAU21SS_vs_OEWI22SS 0.002
## 117 OEAU21SS_vs_OEWI22TW 0.001
## 118 OEAU21SS_vs_OEWI22ML 0.001
## 119 OEAU21SS_vs_OESP22MG 0.001
## 120 OEAU21SS_vs_OESP22BB 0.001
## 121 OEAU21SS_vs_OESP22SS 0.001
## 122 OEAU21SS_vs_OESP22TW 0.001
## 123 OEAU21SS_vs_OESP22ML 0.001
## 124 OEAU21TW_vs_OEAU21ML 0.008
## 125 OEAU21TW_vs_OEWI22BB 0.001
## 126 OEAU21TW_vs_OEWI22MG 0.001
## 127 OEAU21TW_vs_OEWI22SS 0.001
## 128 OEAU21TW_vs_OEWI22TW 0.001
## 129 OEAU21TW_vs_OEWI22ML 0.003
## 130 OEAU21TW_vs_OESP22MG 0.002
## 131 OEAU21TW_vs_OESP22BB 0.002
## 132 OEAU21TW_vs_OESP22SS 0.001
## 133 OEAU21TW_vs_OESP22TW 0.001
## 134 OEAU21TW_vs_OESP22ML 0.001
## 135 OEAU21ML_vs_OEWI22BB 0.002
## 136 OEAU21ML_vs_OEWI22MG 0.002
## 137 OEAU21ML_vs_OEWI22SS 0.001
## 138 OEAU21ML_vs_OEWI22TW 0.003
## 139 OEAU21ML_vs_OEWI22ML 0.002
## 140 OEAU21ML_vs_OESP22MG 0.001
## 141 OEAU21ML_vs_OESP22BB 0.001
## 142 OEAU21ML_vs_OESP22SS 0.001
## 143 OEAU21ML_vs_OESP22TW 0.001
## 144 OEAU21ML_vs_OESP22ML 0.001
## 145 OEWI22BB_vs_OEWI22TW 0.012
## 146 OEWI22BB_vs_OESP22MG 0.001
## 147 OEWI22BB_vs_OESP22BB 0.001
## 148 OEWI22BB_vs_OESP22SS 0.001
## 149 OEWI22BB_vs_OESP22TW 0.001
## 150 OEWI22BB_vs_OESP22ML 0.001
## 151 OEWI22MG_vs_OEWI22TW 0.027
## 152 OEWI22MG_vs_OESP22MG 0.001
## 153 OEWI22MG_vs_OESP22BB 0.001
## 154 OEWI22MG_vs_OESP22SS 0.001
## 155 OEWI22MG_vs_OESP22TW 0.002
## 156 OEWI22MG_vs_OESP22ML 0.002
## 157 OEWI22SS_vs_OEWI22TW 0.045
## 158 OEWI22SS_vs_OESP22MG 0.001
## 159 OEWI22SS_vs_OESP22BB 0.001
## 160 OEWI22SS_vs_OESP22SS 0.001
## 161 OEWI22SS_vs_OESP22TW 0.001
## 162 OEWI22SS_vs_OESP22ML 0.001
## 163 OEWI22TW_vs_OEWI22ML 0.019
## 164 OEWI22TW_vs_OESP22MG 0.001
## 165 OEWI22TW_vs_OESP22BB 0.001
## 166 OEWI22TW_vs_OESP22SS 0.001
## 167 OEWI22TW_vs_OESP22TW 0.003
## 168 OEWI22TW_vs_OESP22ML 0.001
## 169 OEWI22ML_vs_OESP22MG 0.001
## 170 OEWI22ML_vs_OESP22BB 0.001
## 171 OEWI22ML_vs_OESP22SS 0.003
## 172 OEWI22ML_vs_OESP22TW 0.001
## 173 OEWI22ML_vs_OESP22ML 0.001
## 174 OESP22MG_vs_OESP22BB 0.003
## 175 OESP22MG_vs_OESP22SS 0.001
## 176 OESP22MG_vs_OESP22TW 0.002
## 177 OESP22MG_vs_OESP22ML 0.001
## 178 OESP22BB_vs_OESP22SS 0.001
## 179 OESP22BB_vs_OESP22TW 0.001
## 180 OESP22BB_vs_OESP22ML 0.001
## 181 OESP22SS_vs_OESP22TW 0.002
## 182 OESP22SS_vs_OESP22ML 0.039
## 183 OESP22TW_vs_OESP22ML 0.001
## [1] "GC"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates)
## Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 14
## Residual 71
## Total 85
## SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 6.2966
## Residual 10.6271
## Total 16.9237
## R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 0.37206
## Residual 0.62794
## Total 1.00000
## F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 3.0048
## Residual
## Total
## Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 0.001
## Residual
## Total
##
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates ***
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "GC"
## [1] "Trust Adonis"
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 14 0.50195 0.035854 1.5555 0.1141
## Residuals 71 1.63652 0.023050
## [1] "Next analysis"
## Comparison p
## 1 GCSU21BB_vs_GCSU21TW 0.008
## 2 GCSU21BB_vs_GCSU21ML 0.001
## 3 GCSU21BB_vs_GCAU21BB 0.022
## 4 GCSU21BB_vs_GCAU21SS 0.005
## 5 GCSU21BB_vs_GCAU21ML 0.037
## 6 GCSU21BB_vs_GCWI22BB 0.035
## 7 GCSU21BB_vs_GCWI22SS 0.005
## 8 GCSU21BB_vs_GCWI22ML 0.016
## 9 GCSU21BB_vs_GCSP22BB 0.032
## 10 GCSU21BB_vs_GCSP22TW 0.044
## 11 GCSU21BB_vs_GCSP22ML 0.007
## 12 GCSU21SS_vs_GCSU21TW 0.036
## 13 GCSU21SS_vs_GCSU21ML 0.005
## 14 GCSU21SS_vs_GCAU21BB 0.012
## 15 GCSU21SS_vs_GCAU21SS 0.001
## 16 GCSU21SS_vs_GCAU21ML 0.014
## 17 GCSU21SS_vs_GCWI22BB 0.040
## 18 GCSU21SS_vs_GCWI22SS 0.013
## 19 GCSU21SS_vs_GCSP22BB 0.002
## 20 GCSU21SS_vs_GCSP22ML 0.001
## 21 GCSU21TW_vs_GCSU21ML 0.030
## 22 GCSU21TW_vs_GCAU21BB 0.004
## 23 GCSU21TW_vs_GCAU21SS 0.002
## 24 GCSU21TW_vs_GCAU21ML 0.004
## 25 GCSU21TW_vs_GCWI22BB 0.006
## 26 GCSU21TW_vs_GCWI22SS 0.001
## 27 GCSU21TW_vs_GCWI22ML 0.012
## 28 GCSU21TW_vs_GCSP22BB 0.003
## 29 GCSU21TW_vs_GCSP22SS 0.005
## 30 GCSU21TW_vs_GCSP22TW 0.006
## 31 GCSU21TW_vs_GCSP22ML 0.001
## 32 GCSU21ML_vs_GCAU21BB 0.001
## 33 GCSU21ML_vs_GCAU21SS 0.001
## 34 GCSU21ML_vs_GCAU21ML 0.001
## 35 GCSU21ML_vs_GCWI22SS 0.003
## 36 GCSU21ML_vs_GCWI22ML 0.014
## 37 GCSU21ML_vs_GCSP22BB 0.001
## 38 GCSU21ML_vs_GCSP22SS 0.001
## 39 GCSU21ML_vs_GCSP22TW 0.001
## 40 GCSU21ML_vs_GCSP22ML 0.001
## 41 GCAU21BB_vs_GCAU21SS 0.006
## 42 GCAU21BB_vs_GCWI22BB 0.037
## 43 GCAU21BB_vs_GCWI22SS 0.019
## 44 GCAU21BB_vs_GCSP22BB 0.036
## 45 GCAU21BB_vs_GCSP22SS 0.024
## 46 GCAU21BB_vs_GCSP22ML 0.001
## 47 GCAU21SS_vs_GCWI22BB 0.021
## 48 GCAU21SS_vs_GCWI22SS 0.002
## 49 GCAU21SS_vs_GCWI22ML 0.005
## 50 GCAU21SS_vs_GCSP22BB 0.006
## 51 GCAU21SS_vs_GCSP22SS 0.012
## 52 GCAU21SS_vs_GCSP22TW 0.033
## 53 GCAU21SS_vs_GCSP22ML 0.002
## 54 GCAU21TW_vs_GCSP22BB 0.030
## 55 GCAU21TW_vs_GCSP22ML 0.032
## 56 GCAU21ML_vs_GCWI22BB 0.041
## 57 GCAU21ML_vs_GCWI22SS 0.044
## 58 GCAU21ML_vs_GCSP22BB 0.012
## 59 GCWI22BB_vs_GCSP22BB 0.022
## 60 GCWI22BB_vs_GCSP22SS 0.032
## 61 GCWI22SS_vs_GCWI22ML 0.034
## 62 GCWI22SS_vs_GCSP22BB 0.001
## 63 GCWI22SS_vs_GCSP22SS 0.038
## 64 GCWI22SS_vs_GCSP22TW 0.028
## 65 GCWI22ML_vs_GCSP22BB 0.009
## 66 GCWI22ML_vs_GCSP22SS 0.030
## 67 GCWI22ML_vs_GCSP22TW 0.030
## 68 GCWI22ML_vs_GCSP22ML 0.011
## 69 GCSP22BB_vs_GCSP22TW 0.034
## 70 GCSP22BB_vs_GCSP22ML 0.005pslist$PERMANOVA_list <- PERMANOVA_listhttps://astrobiomike.github.io/amplicon/dada2_workflow_ex
Checking by all sample types, we get a significant result (0.002) from the betadisper test. This tells us that there is a difference between group dispersions, which means that we can’t trust the results of an adonis (permutational anova) test on this, because the assumption of homogenous within-group disperions is not met. This isn’t all that surprising considering how different the water and biofilm samples are from the rocks.
An important note: An NMDS is just a visualization technique, and is not a statistical assessment of sample separation or correlation. For this you should run a statistical test, like an ANOSIM test for categorical variables, or a Mantel test for continuous variables. Other ordination techniques like PCA, CA, CCA, RDA, etc, may be more useful to you if you want your axes to be meaningful, or if you want to talk about variation partitioning.
6.3.4 dbRDA / CAP
distance based RDA automatically removing Oxygen values for multicollinearity with almost all other abiotic measurements, further there is an automatic exclusion of other multicollinear factors, inspect the output datasets instead of trusting the RDA plot from scratch
RDA_list <- list()
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
require(phyloseq)
require(vegan)
require(plyr)
require(dplyr)
require(ggrepel)
require(cowplot)
require(ggordiplots)
Datasetname <- sub("ps_", "", Dataset)
ps <- pslist[[Dataset]]
RDA_list_length <- length(RDA_list)
GroupingVariables <- c("fSeason", "fLOC")
if (Datasetname %in% c("GC", "OE")) {
Traits <- c(
"Salinity", "Temperature",
"Age", "FCF", "HSI", "SSI", "GSI", "FillLevel",
"NH4", "NO2", "NO3", "O2", "PO4", "pH", "SPM", "TOC",
GroupingVariables)
print(Traits)
} else if (Datasetname %in% c("WF")) {
Traits <- c(
"Salinity", "Temperature", "NH4", "NO2", "NO3", "O2", "PO4", "pH", "SPM", "TOC", GroupingVariables)
print(Traits)
} else {print("error")}
if (Datasetname %in% c("GC", "OE")) {
print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
} else if (Datasetname %in% c("WF")) {
print(paste("No Waterfiler Samples removed ;", Dataset, sep = ""))
ps_Filter <- ps
}
'%!in%' <- function(x,y)!('%in%'(x,y))
Include <- Traits[Traits %!in% GroupingVariables]
df_Traits <- data.frame(sample_data(ps_Filter))
df_Traits <- df_Traits[names(df_Traits) %in% Traits]
df_Include <- df_Traits[names(df_Traits) %in% Include]
print("NAs in TraitData")
print(table(is.na(df_Include)))
print((df_Include[!complete.cases(df_Include), ]))
df_Include <- df_Include[complete.cases(df_Include), ]
df_Traits <- df_Traits[complete.cases(df_Traits), ]
dim(df_Traits)
df_Include <- decostand(df_Include, method='standardize')
dim(df_Include)
ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
rownames(df_Include[complete.cases(df_Include), ])], ps_Filter)
##################################################
#Determine which variables are useful for the plot#
##################################################
abund_table <- as.data.frame(otu_table(ps_Filter, taxa_are_rows = FALSE))
print("Smallest Group for avgdist")
smallest_group <- min(rowSums(otu_table(ps_Filter)))
print(smallest_group)
avgdist_abund_table <- avgdist(abund_table, dmethod = "bray", sample = smallest_group)
dpRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Include, distance='bray')
abund_table.adonis <- vegan::adonis2(avgdist_abund_table ~ ., data=df_Include)
print(paste("PERMANOVA", Datasetname))
print(abund_table.adonis)
abund_table.adonis.sig <- abund_table.adonis[!is.na(abund_table.adonis$`Pr(>F)`),]
abund_table.adonis.sig <- abund_table.adonis.sig[abund_table.adonis.sig$`Pr(>F)` < 0.05,]
PERMANOVAsig.variables <- row.names(abund_table.adonis.sig)
###########################################################################
#Compute Distance-based Redundancy Analysis based on average bray distance#
###########################################################################
dbRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Include, distance='bray')
dbrda_plot = plot(dpRDA_results, xlab=ord_labels(dpRDA_results)[1], ylab=ord_labels(dpRDA_results)[2])
anova(dbRDA_results) # permutation test of statistical significance
####################
#Variable selection#
####################
######################
#Check for redundancy#
######################
print("VIF")
print(sort(vif.cca(dpRDA_results))) # calculate VIF, scores >10 are redundant
vars.envfit <- names(which(vif.cca(dpRDA_results) <= 10))
##################
# model selection#
##################
rda1 <- capscale(avgdist_abund_table ~ ., data=df_Include) # set up full and null models for 'ordistep' full model
rda0 <- capscale(avgdist_abund_table ~ 1, data=df_Include) # intercept-only (null) model
step.env <- ordistep(rda0, scope=formula(rda1), direction='both')
#### # perform forward and backward selection of explanatory variables
print(paste("Ordisep", Datasetname))
print(step.env$anova)
vars.ordistep <- gsub('^. ', '', rownames(step.env$anova)) # get variable names from 'ordistep' and 'envfit' results
Significant <- unique(vars.ordistep)
df_Significant <- df_Traits[names(df_Traits) %in% vars.ordistep]
#First check of O2 is collinear, otherwise VIF goes up for many values
if ("O2" %in% names(which(vif.cca(dbRDA_results) > 10))) {
print(paste("Removing O2 for Collinearity"))
df_Significant <- df_Significant[, !names(df_Significant) %in% c("O2")]
}
# Perform the same RDA but with fewer variables (selected/relevant variables)
dbRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Significant, distance='bray')
summary(dbRDA_results, display=NULL) # summary of the results
anova(dbRDA_results) # permutation test of statistical significance
print("VIF")
print(sort(vif.cca(dbRDA_results))) # calculate VIF, scores >10 are redundant
vars.envfit <- names(which(vif.cca(dbRDA_results) <= 10))
if (length(names(which(vif.cca(dbRDA_results) > 10))) > 0) {
print(paste("Removing", names(which(vif.cca(dbRDA_results) > 10)), "for Collinearity"))
df_Significant <- df_Significant[, !names(df_Significant) %in% names(which(vif.cca(dbRDA_results) > 10))]
}
# Perform the same RDA but with fewer variables (selected/relevant variables)
dbRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Significant, distance='bray')
summary(dbRDA_results, display=NULL) # summary of the results
anova(dbRDA_results) # permutation test of statistical significance
envfit_res <- envfit(dbRDA_results, df_Significant, permutations = 999, na.rm = TRUE)
RDA_list[[RDA_list_length+1]] <- dbRDA_results
names(RDA_list)[[RDA_list_length+1]] <- paste("dbRDA_results", Datasetname, sep="_")
RDA_list[[RDA_list_length+2]] <- step.env
names(RDA_list)[[RDA_list_length+2]] <- paste("ordistep_results", Datasetname, sep="_")
RDA_list[[RDA_list_length+3]] <- abund_table.adonis.sig
names(RDA_list)[[RDA_list_length+3]] <- paste("PERMANOVA_results", Datasetname, sep="_")
RDA_list[[RDA_list_length+4]] <- envfit_res
names(RDA_list)[[RDA_list_length+4]] <- paste("EnvFit_results", Datasetname, sep="_")
# Set up dataframes for plotting the results
scores_df <- cbind(df_Traits, scores(dbRDA_results, display="sites"))
#spp_dpRDA <- cbind(data.frame(tax_table(ps_Filter)), vegan::scores(dbRDA_results, display="species"))
scores_dbRDA <- vegan::scores(dbRDA_results, display=c("sp","wa","lc","bp","cn"))
##########
#Plot RDA#
##########
tryCatch({
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% sample_data(ps)$Season]
p <-ggplot(scores_df, aes(x=CAP1, y=CAP2, color=fSeason, shape=fLOC)) + #, shape=RepsSpecs
ggtitle (paste(Datasetname, "RDA"))+
scale_color_manual(values = COL) +
scale_fill_manual(values = COL) +
geom_point(size=4.5)
if (Datasetname %in% c("GC")) {
p <- p + scale_shape_manual(values=c(16, 17, 8, 25))
} else if (Datasetname %in% c("OE", "WF")) {
p <- p + scale_shape_manual(values=c(15, 16, 17, 8, 25))
}
#p <- p + ggrepel::geom_text_repel(aes(label = fLOC))
multiplier <- vegan:::ordiArrowMul(scores_dbRDA$biplot)
df_arrows <- scores_dbRDA$biplot*multiplier
colnames(df_arrows)<-c("x","y")
df_arrows=as.data.frame(df_arrows)
p <- p + geom_segment(data=df_arrows,inherit.aes=F, mapping=(aes(x = 0, y = 0, xend = x, yend = y)),
arrow = arrow(length = unit(0.2, "cm")),color="grey60",alpha=0.8) +
geom_text_repel(data=as.data.frame(df_arrows*1.2), inherit.aes=F, mapping=aes(x, y, label =
rownames(df_arrows)),color="grey20",alpha=0.99, size = 6)
p <- p + atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=10),
legend.text=element_text(size=10)) +
theme(
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "white"))
# add in Percent to Axis lables
p <- p + xlab(paste("RDA1", sub("CAP1", "", ord_labels(dbRDA_results)[1]))) +
ylab(paste("RDA2", sub("CAP2", "", ord_labels(dbRDA_results)[2]))) +
theme(
axis.title.x.bottom = element_text(size=12,face = "bold"),
axis.title.y.left = element_text(size=12,face = "bold"))
# p <- p + annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1, size = 5,
# label = paste(
# paste(
# "p:", pslist$PERMANOVA_list[[Dataset]]$Season$Adonis$`Pr(>F)`[1]),
#
# paste(
# "DF:", pslist$PERMANOVA_list[[Dataset]]$Season$Adonis$Df[3]),
#
# paste(
# "F:",round(pslist$PERMANOVA_list[[Dataset]]$Season$Adonis$F[1])),
# sep = "\n"),
# color = c("darkred", "darkred", "darkred"), fontface = "bold")
RDA_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
ggsave(RDA_plot, filename = paste(paste(save_name, "dbRDA", Datasetname, sep="_"), ".png", sep=""),
path = pathPlots, device='png', dpi=300, width = 7,height = 5)
},error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
plot(RDA_plot)
RDA_list[[RDA_list_length+5]] <- RDA_plot
names(RDA_list)[[RDA_list_length+5]] <- paste("dbRDA", Datasetname, sep="_")
}## [1] "Salinity" "Temperature" "Age" "FCF" "HSI"
## [6] "SSI" "GSI" "FillLevel" "NH4" "NO2"
## [11] "NO3" "O2" "PO4" "pH" "SPM"
## [16] "TOC" "fSeason" "fLOC"
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OE"
## named integer(0)
## [1] "NAs in TraitData"
##
## FALSE
## 2208
## [1] Temperature Salinity FillLevel Age HSI SSI
## [7] GSI FCF NH4 NO2 NO3 O2
## [13] PO4 pH TOC SPM
## <0 rows> (or 0-length row.names)
## [1] "Smallest Group for avgdist"
## [1] 9071
## [1] "PERMANOVA OE"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = avgdist_abund_table ~ ., data = df_Include)
## Df SumOfSqs R2 F Pr(>F)
## Temperature 1 3.333 0.10234 25.3357 0.001 ***
## Salinity 1 1.123 0.03448 8.5365 0.001 ***
## FillLevel 1 0.550 0.01689 4.1807 0.002 **
## Age 1 0.290 0.00892 2.2070 0.032 *
## HSI 1 1.185 0.03639 9.0099 0.001 ***
## SSI 1 0.361 0.01108 2.7428 0.013 *
## GSI 1 3.260 0.10009 24.7785 0.001 ***
## FCF 1 0.597 0.01832 4.5345 0.001 ***
## NH4 1 0.477 0.01463 3.6229 0.007 **
## NO2 1 0.911 0.02797 6.9252 0.001 ***
## NO3 1 0.910 0.02794 6.9165 0.001 ***
## O2 1 0.695 0.02133 5.2796 0.001 ***
## PO4 1 1.146 0.03519 8.7124 0.001 ***
## pH 1 0.558 0.01713 4.2398 0.002 **
## TOC 1 0.667 0.02047 5.0665 0.001 ***
## SPM 1 0.588 0.01807 4.4725 0.001 ***
## Residual 121 15.920 0.48877
## Total 137 32.571 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "VIF"
## SSI FillLevel Age SPM HSI PO4
## 1.235635 1.238250 1.389237 1.907918 2.037286 2.619247
## FCF TOC GSI pH Salinity NH4
## 2.658719 2.989917 3.585188 5.062142 5.551753 7.980814
## NO2 NO3 Temperature O2
## 8.216231 8.224014 22.591065 27.194257
##
## Start: avgdist_abund_table ~ 1
##
## Df AIC F Pr(>F)
## + GSI 1 460.78 28.1692 0.005 **
## + Temperature 1 472.18 15.1623 0.005 **
## + FCF 1 474.40 12.7475 0.005 **
## + SPM 1 475.32 11.7572 0.005 **
## + O2 1 477.93 8.9884 0.005 **
## + NO3 1 479.03 7.8392 0.005 **
## + TOC 1 479.25 7.6061 0.005 **
## + PO4 1 479.68 7.1594 0.005 **
## + HSI 1 479.71 7.1281 0.005 **
## + Salinity 1 481.71 5.0705 0.005 **
## + pH 1 482.40 4.3649 0.005 **
## + NO2 1 483.21 3.5450 0.005 **
## + Age 1 484.02 2.7252 0.010 **
## + NH4 1 484.46 2.2894 0.055 .
## + FillLevel 1 484.76 1.9892 0.055 .
## + SSI 1 485.50 1.2535 0.200
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI
##
## Df AIC F Pr(>F)
## - GSI 1 484.76 28.169 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Temperature 1 452.19 10.7722 0.005 **
## + O2 1 453.80 9.0860 0.005 **
## + NO3 1 454.29 8.5733 0.005 **
## + FCF 1 456.09 6.7056 0.005 **
## + TOC 1 456.51 6.2787 0.005 **
## + Salinity 1 457.14 5.6404 0.005 **
## + PO4 1 457.21 5.5660 0.005 **
## + SPM 1 457.73 5.0328 0.005 **
## + NO2 1 458.00 4.7602 0.005 **
## + pH 1 458.67 4.0854 0.005 **
## + NH4 1 459.92 2.8358 0.005 **
## + FillLevel 1 459.76 2.9901 0.015 *
## + HSI 1 460.81 1.9472 0.035 *
## + Age 1 461.18 1.5752 0.080 .
## + SSI 1 461.85 0.9143 0.455
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature
##
## Df AIC F Pr(>F)
## - Temperature 1 460.78 10.772 0.005 **
## - GSI 1 472.18 23.315 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + O2 1 445.64 8.5701 0.005 **
## + PO4 1 448.07 6.0803 0.005 **
## + Salinity 1 448.14 6.0099 0.005 **
## + SPM 1 448.32 5.8262 0.005 **
## + NO2 1 449.13 5.0010 0.005 **
## + pH 1 449.59 4.5392 0.005 **
## + TOC 1 450.14 3.9902 0.005 **
## + NO3 1 450.19 3.9417 0.005 **
## + FCF 1 450.70 3.4284 0.005 **
## + NH4 1 451.03 3.1051 0.005 **
## + FillLevel 1 451.07 3.0616 0.005 **
## + HSI 1 451.86 2.2801 0.015 *
## + Age 1 452.71 1.4465 0.165
## + SSI 1 453.72 0.4526 0.915
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2
##
## Df AIC F Pr(>F)
## - O2 1 452.19 8.5701 0.005 **
## - Temperature 1 453.80 10.2385 0.005 **
## - GSI 1 463.71 20.9806 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + pH 1 439.44 8.1328 0.005 **
## + Salinity 1 441.28 6.2722 0.005 **
## + SPM 1 441.30 6.2451 0.005 **
## + TOC 1 442.83 4.7153 0.005 **
## + NO3 1 443.28 4.2677 0.005 **
## + FCF 1 443.70 3.8510 0.005 **
## + PO4 1 443.91 3.6382 0.005 **
## + NO2 1 444.71 2.8483 0.010 **
## + FillLevel 1 445.02 2.5461 0.010 **
## + NH4 1 445.28 2.2890 0.010 **
## + Age 1 446.19 1.4033 0.140
## + HSI 1 446.50 1.1024 0.290
## + SSI 1 447.03 0.5882 0.780
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH
##
## Df AIC F Pr(>F)
## - pH 1 445.64 8.1328 0.005 **
## - O2 1 449.59 12.2392 0.005 **
## - Temperature 1 450.23 12.9127 0.005 **
## - GSI 1 453.85 16.7908 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + SPM 1 434.74 6.5718 0.005 **
## + Salinity 1 436.49 4.8270 0.005 **
## + PO4 1 437.22 4.1023 0.005 **
## + FCF 1 438.11 3.2260 0.005 **
## + NO3 1 438.45 2.8928 0.005 **
## + NH4 1 438.75 2.6028 0.005 **
## + NO2 1 438.76 2.5914 0.005 **
## + TOC 1 438.76 2.5889 0.010 **
## + FillLevel 1 439.57 1.8039 0.035 *
## + Age 1 439.92 1.4699 0.085 .
## + HSI 1 440.31 1.0943 0.290
## + SSI 1 440.77 0.6427 0.800
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM
##
## Df AIC F Pr(>F)
## - SPM 1 439.44 6.5718 0.005 **
## - pH 1 441.30 8.4503 0.005 **
## - O2 1 445.43 12.7130 0.005 **
## - Temperature 1 446.34 13.6720 0.005 **
## - GSI 1 447.19 14.5702 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + PO4 1 432.19 4.3880 0.005 **
## + Salinity 1 432.31 4.2773 0.005 **
## + NO3 1 433.65 2.9648 0.005 **
## + FCF 1 434.10 2.5304 0.005 **
## + TOC 1 433.78 2.8350 0.010 **
## + NO2 1 434.94 1.7204 0.035 *
## + NH4 1 434.84 1.8165 0.045 *
## + FillLevel 1 434.94 1.7228 0.065 .
## + HSI 1 435.56 1.1214 0.235
## + Age 1 435.52 1.1629 0.260
## + SSI 1 436.04 0.6694 0.750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4
##
## Df AIC F Pr(>F)
## - PO4 1 434.74 4.3880 0.010 **
## - SPM 1 437.22 6.8446 0.005 **
## - pH 1 439.27 8.9097 0.005 **
## - O2 1 441.19 10.8699 0.005 **
## - GSI 1 442.52 12.2374 0.005 **
## - Temperature 1 443.00 12.7440 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + NO3 1 428.65 5.3231 0.005 **
## + Salinity 1 430.59 3.4383 0.005 **
## + TOC 1 430.64 3.3855 0.005 **
## + FCF 1 431.53 2.5314 0.005 **
## + NO2 1 432.35 1.7433 0.020 *
## + FillLevel 1 432.38 1.7159 0.035 *
## + NH4 1 432.58 1.5329 0.065 .
## + Age 1 432.98 1.1488 0.215
## + HSI 1 433.08 1.0523 0.325
## + SSI 1 433.56 0.5969 0.835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3
##
## Df AIC F Pr(>F)
## - NO3 1 432.19 5.3231 0.005 **
## - pH 1 432.95 6.0722 0.005 **
## - PO4 1 433.65 6.7608 0.005 **
## - SPM 1 434.12 7.2311 0.005 **
## - GSI 1 434.56 7.6664 0.005 **
## - O2 1 435.75 8.8595 0.005 **
## - Temperature 1 440.64 13.8615 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Salinity 1 423.44 6.9246 0.005 **
## + FCF 1 428.22 2.2972 0.005 **
## + TOC 1 428.13 2.3805 0.020 *
## + NH4 1 428.98 1.5759 0.040 *
## + NO2 1 428.79 1.7537 0.045 *
## + FillLevel 1 429.03 1.5267 0.090 .
## + Age 1 429.38 1.1964 0.225
## + HSI 1 429.49 1.0972 0.360
## + SSI 1 429.77 0.8314 0.640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity
##
## Df AIC F Pr(>F)
## - GSI 1 426.94 5.2419 0.005 **
## - SPM 1 427.51 5.7997 0.005 **
## - pH 1 427.78 6.0653 0.005 **
## - Salinity 1 428.65 6.9246 0.005 **
## - PO4 1 429.00 7.2697 0.005 **
## - NO3 1 430.59 8.8445 0.005 **
## - O2 1 431.10 9.3589 0.005 **
## - Temperature 1 437.01 15.4045 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + TOC 1 419.53 5.5965 0.005 **
## + NH4 1 423.47 1.8371 0.015 *
## + NO2 1 423.42 1.8848 0.020 *
## + FCF 1 423.51 1.8049 0.030 *
## + FillLevel 1 423.77 1.5582 0.120
## + Age 1 424.13 1.2183 0.210
## + HSI 1 424.21 1.1473 0.310
## + SSI 1 424.47 0.8999 0.575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC
##
## Df AIC F Pr(>F)
## - GSI 1 420.82 3.0869 0.020 *
## - TOC 1 423.44 5.5965 0.005 **
## - SPM 1 423.74 5.8872 0.005 **
## - PO4 1 424.11 6.2432 0.005 **
## - pH 1 425.51 7.6131 0.005 **
## - NO3 1 425.81 7.9113 0.005 **
## - O2 1 427.14 9.2248 0.005 **
## - Salinity 1 428.13 10.2172 0.005 **
## - Temperature 1 434.06 16.2840 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + NH4 1 419.22 2.1473 0.015 *
## + NO2 1 419.43 1.9469 0.015 *
## + FillLevel 1 419.83 1.5754 0.060 .
## + Age 1 420.18 1.2535 0.145
## + HSI 1 420.24 1.1922 0.215
## + FCF 1 420.41 1.0369 0.325
## + SSI 1 420.73 0.7391 0.670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4
##
## Df AIC F Pr(>F)
## - NH4 1 419.53 2.1473 0.055 .
## - GSI 1 420.17 2.7486 0.015 *
## - PO4 1 422.01 4.4859 0.005 **
## - SPM 1 422.29 4.7550 0.005 **
## - TOC 1 423.47 5.8868 0.005 **
## - NO3 1 425.23 7.5899 0.005 **
## - O2 1 425.42 7.7784 0.005 **
## - pH 1 425.57 7.9256 0.005 **
## - Salinity 1 428.08 10.4009 0.005 **
## - Temperature 1 430.99 13.3241 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + NO2 1 417.65 3.3064 0.005 **
## + FillLevel 1 419.49 1.5879 0.030 *
## + Age 1 419.85 1.2598 0.170
## + HSI 1 419.91 1.1990 0.205
## + FCF 1 420.20 0.9310 0.560
## + SSI 1 420.41 0.7395 0.685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2
##
## Df AIC F Pr(>F)
## - GSI 1 418.68 2.7972 0.020 *
## - NO2 1 419.22 3.3064 0.005 **
## - NH4 1 419.43 3.5074 0.005 **
## - PO4 1 420.70 4.7040 0.005 **
## - SPM 1 420.83 4.8218 0.005 **
## - O2 1 421.78 5.7240 0.005 **
## - TOC 1 422.03 5.9651 0.005 **
## - pH 1 422.73 6.6352 0.005 **
## - NO3 1 423.90 7.7626 0.005 **
## - Salinity 1 426.74 10.5452 0.005 **
## - Temperature 1 428.46 12.2570 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + FillLevel 1 417.62 1.8441 0.020 *
## + HSI 1 418.07 1.4305 0.100 .
## + Age 1 418.29 1.2315 0.205
## + FCF 1 418.60 0.9538 0.440
## + SSI 1 418.71 0.8502 0.590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel
##
## Df AIC F Pr(>F)
## - FillLevel 1 417.65 1.8441 0.060 .
## - GSI 1 418.70 2.8183 0.015 *
## - NO2 1 419.49 3.5527 0.005 **
## - NH4 1 419.72 3.7627 0.005 **
## - SPM 1 420.46 4.4568 0.005 **
## - PO4 1 420.68 4.6652 0.005 **
## - O2 1 421.09 5.0512 0.005 **
## - TOC 1 422.03 5.9418 0.005 **
## - pH 1 422.11 6.0116 0.005 **
## - NO3 1 423.79 7.6171 0.005 **
## - Salinity 1 426.73 10.4791 0.005 **
## - Temperature 1 427.44 11.1726 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + HSI 1 417.91 1.5464 0.050 *
## + Age 1 418.19 1.2916 0.140
## + FCF 1 418.61 0.9150 0.550
## + SSI 1 418.75 0.7909 0.645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel + HSI
##
## Df AIC F Pr(>F)
## - HSI 1 417.62 1.5464 0.125
## - FillLevel 1 418.07 1.9570 0.055 .
## - GSI 1 418.72 2.5434 0.015 *
## - NH4 1 420.22 3.9314 0.010 **
## - NO2 1 420.16 3.8743 0.005 **
## - O2 1 420.77 4.4393 0.005 **
## - SPM 1 420.77 4.4451 0.005 **
## - PO4 1 421.02 4.6751 0.005 **
## - TOC 1 422.37 5.9367 0.005 **
## - pH 1 422.45 6.0105 0.005 **
## - NO3 1 424.24 7.7127 0.005 **
## - Temperature 1 426.68 10.0653 0.005 **
## - Salinity 1 427.16 10.5260 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel
##
## Df AIC F Pr(>F)
## + HSI 1 417.91 1.5464 0.045 *
## + Age 1 418.19 1.2916 0.215
## + FCF 1 418.61 0.9150 0.525
## + SSI 1 418.75 0.7909 0.605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel + HSI
##
## Df AIC F Pr(>F)
## - HSI 1 417.62 1.5464 0.105
## - FillLevel 1 418.07 1.9570 0.065 .
## - GSI 1 418.72 2.5434 0.025 *
## - NO2 1 420.16 3.8743 0.010 **
## - NH4 1 420.22 3.9314 0.005 **
## - O2 1 420.77 4.4393 0.005 **
## - SPM 1 420.77 4.4451 0.005 **
## - PO4 1 421.02 4.6751 0.005 **
## - TOC 1 422.37 5.9367 0.005 **
## - pH 1 422.45 6.0105 0.005 **
## - NO3 1 424.24 7.7127 0.005 **
## - Temperature 1 426.68 10.0653 0.005 **
## - Salinity 1 427.16 10.5260 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel
##
## Df AIC F Pr(>F)
## + HSI 1 417.91 1.5464 0.050 *
## + Age 1 418.19 1.2916 0.125
## + FCF 1 418.61 0.9150 0.505
## + SSI 1 418.75 0.7909 0.690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel + HSI
##
## Df AIC F Pr(>F)
## - HSI 1 417.62 1.5464 0.145
## - FillLevel 1 418.07 1.9570 0.075 .
## - GSI 1 418.72 2.5434 0.035 *
## - NO2 1 420.16 3.8743 0.005 **
## - NH4 1 420.22 3.9314 0.005 **
## - O2 1 420.77 4.4393 0.005 **
## - SPM 1 420.77 4.4451 0.005 **
## - PO4 1 421.02 4.6751 0.005 **
## - TOC 1 422.37 5.9367 0.005 **
## - pH 1 422.45 6.0105 0.005 **
## - NO3 1 424.24 7.7127 0.005 **
## - Temperature 1 426.68 10.0653 0.005 **
## - Salinity 1 427.16 10.5260 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 + NO3 + Salinity + TOC + NH4 + NO2 + FillLevel
##
## Df AIC F Pr(>F)
## + HSI 1 417.91 1.5464 0.060 .
## + Age 1 418.19 1.2916 0.125
## + FCF 1 418.61 0.9150 0.525
## + SSI 1 418.75 0.7909 0.680
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## - FillLevel 1 417.65 1.8441 0.075 .
## - GSI 1 418.70 2.8183 0.010 **
## - NO2 1 419.49 3.5527 0.005 **
## - NH4 1 419.72 3.7627 0.005 **
## - SPM 1 420.46 4.4568 0.005 **
## - PO4 1 420.68 4.6652 0.005 **
## - O2 1 421.09 5.0512 0.005 **
## - TOC 1 422.03 5.9418 0.005 **
## - pH 1 422.11 6.0116 0.005 **
## - NO3 1 423.79 7.6171 0.005 **
## - Salinity 1 426.73 10.4791 0.005 **
## - Temperature 1 427.44 11.1726 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + HSI 1 417.91 1.5464 0.060 .
## + Age 1 418.19 1.2916 0.145
## + FCF 1 418.61 0.9150 0.580
## + SSI 1 418.75 0.7909 0.640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] "Ordisep OE"
## Df AIC F Pr(>F)
## + GSI 1 460.78 28.1692 0.005 **
## + Temperature 1 452.19 10.7722 0.005 **
## + O2 1 445.64 8.5701 0.005 **
## + pH 1 439.44 8.1328 0.005 **
## + SPM 1 434.74 6.5718 0.005 **
## + PO4 1 432.19 4.3880 0.005 **
## + NO3 1 428.65 5.3231 0.005 **
## + Salinity 1 423.44 6.9246 0.005 **
## + TOC 1 419.53 5.5965 0.005 **
## + NH4 1 419.22 2.1473 0.015 *
## + NO2 1 417.65 3.3064 0.005 **
## + FillLevel 1 417.62 1.8441 0.020 *
## + HSI 1 417.91 1.5464 0.050 *
## - HSI 1 417.62 1.5464 0.125
## + HSI1 1 417.91 1.5464 0.045 *
## - HSI1 1 417.62 1.5464 0.105
## + HSI2 1 417.91 1.5464 0.050 *
## - HSI2 1 417.62 1.5464 0.145
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Removing O2 for Collinearity"
## [1] "VIF"
## FillLevel SPM HSI PO4 TOC GSI
## 1.135031 1.471549 1.778025 2.008939 2.489322 2.839159
## pH Salinity Temperature NO2 NO3 NH4
## 3.096458 5.051252 6.326722 6.891286 6.903946 7.898809## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## [1] "Salinity" "Temperature" "Age" "FCF" "HSI"
## [6] "SSI" "GSI" "FillLevel" "NH4" "NO2"
## [11] "NO3" "O2" "PO4" "pH" "SPM"
## [16] "TOC" "fSeason" "fLOC"
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GC"
## GCAU21TWFL2 GCAU21MLEB1
## 45 46
## [1] "NAs in TraitData"
##
## FALSE
## 1376
## [1] Temperature Salinity FillLevel Age HSI SSI
## [7] GSI FCF NH4 NO2 NO3 O2
## [13] PO4 pH TOC SPM
## <0 rows> (or 0-length row.names)
## [1] "Smallest Group for avgdist"
## [1] 5998
## [1] "PERMANOVA GC"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = avgdist_abund_table ~ ., data = df_Include)
## Df SumOfSqs R2 F Pr(>F)
## Temperature 1 0.7239 0.04278 4.6766 0.001 ***
## Salinity 1 0.7376 0.04359 4.7650 0.001 ***
## FillLevel 1 0.1426 0.00843 0.9214 0.491
## Age 1 0.3107 0.01836 2.0070 0.026 *
## HSI 1 0.1643 0.00971 1.0615 0.339
## SSI 1 0.0852 0.00504 0.5504 0.908
## GSI 1 0.1150 0.00680 0.7429 0.719
## FCF 1 0.1592 0.00941 1.0287 0.377
## NH4 1 0.5826 0.03444 3.7639 0.001 ***
## NO2 1 1.0029 0.05927 6.4789 0.001 ***
## NO3 1 0.5381 0.03181 3.4766 0.001 ***
## O2 1 0.2650 0.01566 1.7118 0.052 .
## PO4 1 0.5173 0.03057 3.3418 0.001 ***
## pH 1 0.1817 0.01074 1.1737 0.264
## TOC 1 0.4240 0.02506 2.7392 0.005 **
## SPM 1 0.2886 0.01706 1.8646 0.045 *
## Residual 69 10.6804 0.63127
## Total 85 16.9190 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "VIF"
## SSI HSI FillLevel FCF GSI Age
## 1.333285 1.637703 1.972504 2.141985 2.198489 2.573010
## TOC SPM PO4 Salinity NO3 pH
## 2.964841 3.298999 5.820725 6.755114 6.913369 7.086578
## NO2 NH4 Temperature O2
## 13.252686 18.680343 43.279216 48.706396
##
## Start: avgdist_abund_table ~ 1
##
## Df AIC F Pr(>F)
## + NO2 1 240.77 6.1761 0.005 **
## + PO4 1 241.76 5.1427 0.005 **
## + O2 1 242.06 4.8424 0.005 **
## + NH4 1 242.83 4.0428 0.005 **
## + Temperature 1 243.14 3.7277 0.005 **
## + Salinity 1 243.25 3.6144 0.005 **
## + SPM 1 243.32 3.5498 0.005 **
## + TOC 1 244.10 2.7520 0.010 **
## + Age 1 244.36 2.4953 0.010 **
## + pH 1 244.09 2.7618 0.015 *
## + NO3 1 244.48 2.3697 0.015 *
## + GSI 1 244.95 1.9058 0.055 .
## + HSI 1 245.71 1.1420 0.290
## + FillLevel 1 246.01 0.8504 0.575
## + SSI 1 245.96 0.8968 0.580
## + FCF 1 246.12 0.7377 0.720
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2
##
## Df AIC F Pr(>F)
## - NO2 1 244.88 6.1761 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + PO4 1 236.46 6.3245 0.005 **
## + O2 1 239.22 3.5021 0.005 **
## + Temperature 1 239.22 3.5001 0.005 **
## + TOC 1 239.50 3.2160 0.005 **
## + NO3 1 240.00 2.7179 0.005 **
## + SPM 1 240.06 2.6559 0.005 **
## + Salinity 1 240.26 2.4608 0.015 *
## + pH 1 240.25 2.4675 0.020 *
## + NH4 1 240.35 2.3687 0.020 *
## + Age 1 240.50 2.2198 0.020 *
## + GSI 1 241.10 1.6293 0.095 .
## + HSI 1 241.53 1.2127 0.220
## + FillLevel 1 241.70 1.0425 0.415
## + FCF 1 242.01 0.7424 0.700
## + SSI 1 242.04 0.7132 0.705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2 + PO4
##
## Df AIC F Pr(>F)
## - PO4 1 240.77 6.3245 0.005 **
## - NO2 1 241.76 7.3601 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + TOC 1 235.16 3.2023 0.005 **
## + Temperature 1 235.23 3.1358 0.005 **
## + pH 1 235.55 2.8241 0.010 **
## + O2 1 235.62 2.7539 0.010 **
## + NO3 1 235.64 2.7311 0.010 **
## + Salinity 1 235.86 2.5150 0.015 *
## + Age 1 235.99 2.3840 0.015 *
## + SPM 1 236.24 2.1419 0.015 *
## + GSI 1 236.78 1.6143 0.035 *
## + NH4 1 236.84 1.5590 0.085 .
## + FillLevel 1 237.10 1.3021 0.120
## + HSI 1 237.22 1.1887 0.260
## + SSI 1 237.72 0.7051 0.715
## + FCF 1 237.71 0.7135 0.810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC
##
## Df AIC F Pr(>F)
## - TOC 1 236.46 3.2023 0.005 **
## - PO4 1 239.50 6.2742 0.005 **
## - NO2 1 241.71 8.5623 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Salinity 1 233.51 3.5185 0.005 **
## + pH 1 233.68 3.3516 0.005 **
## + O2 1 234.22 2.8225 0.005 **
## + Age 1 234.55 2.5006 0.005 **
## + NO3 1 234.70 2.3561 0.005 **
## + Temperature 1 233.85 3.1855 0.010 **
## + SPM 1 235.12 1.9481 0.035 *
## + GSI 1 235.37 1.7065 0.075 .
## + HSI 1 235.87 1.2312 0.215
## + NH4 1 235.84 1.2544 0.235
## + FillLevel 1 235.94 1.1571 0.295
## + FCF 1 236.37 0.7551 0.730
## + SSI 1 236.52 0.6071 0.850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity
##
## Df AIC F Pr(>F)
## - TOC 1 235.86 4.2059 0.010 **
## - Salinity 1 235.16 3.5185 0.005 **
## - PO4 1 238.44 6.8023 0.005 **
## - NO2 1 238.46 6.8193 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Temperature 1 231.72 3.6029 0.005 **
## + O2 1 232.39 2.9503 0.005 **
## + pH 1 232.59 2.7648 0.005 **
## + NO3 1 233.07 2.3034 0.010 **
## + SPM 1 233.43 1.9550 0.010 **
## + Age 1 233.03 2.3350 0.015 *
## + GSI 1 233.57 1.8222 0.060 .
## + NH4 1 234.13 1.2891 0.135
## + HSI 1 234.12 1.2965 0.175
## + FillLevel 1 234.15 1.2734 0.210
## + FCF 1 234.77 0.6894 0.860
## + SSI 1 234.86 0.6041 0.865
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity + Temperature
##
## Df AIC F Pr(>F)
## - Temperature 1 233.51 3.6029 0.005 **
## - Salinity 1 233.85 3.9336 0.005 **
## - TOC 1 234.41 4.4878 0.005 **
## - NO2 1 236.18 6.2430 0.005 **
## - PO4 1 236.52 6.5854 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + NO3 1 230.67 2.8494 0.010 **
## + Age 1 231.25 2.2975 0.015 *
## + NH4 1 231.79 1.7909 0.025 *
## + O2 1 231.94 1.6460 0.055 .
## + SPM 1 232.26 1.3523 0.160
## + pH 1 232.47 1.1578 0.260
## + FillLevel 1 232.61 1.0224 0.380
## + HSI 1 232.62 1.0152 0.445
## + FCF 1 232.85 0.7977 0.695
## + GSI 1 233.00 0.6620 0.820
## + SSI 1 233.06 0.6115 0.900
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity + Temperature + NO3
##
## Df AIC F Pr(>F)
## - NO3 1 231.72 2.8494 0.005 **
## - TOC 1 232.48 3.5741 0.005 **
## - Salinity 1 232.69 3.7767 0.005 **
## - PO4 1 232.87 3.9555 0.005 **
## - Temperature 1 233.07 4.1417 0.005 **
## - NO2 1 234.53 5.5672 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Age 1 230.67 1.8337 0.035 *
## + SPM 1 231.17 1.3753 0.080 .
## + NH4 1 230.97 1.5560 0.115
## + O2 1 231.19 1.3597 0.170
## + pH 1 231.41 1.1531 0.275
## + HSI 1 231.48 1.0886 0.365
## + FillLevel 1 231.71 0.8739 0.575
## + FCF 1 231.84 0.7592 0.635
## + GSI 1 231.93 0.6804 0.855
## + SSI 1 232.03 0.5823 0.930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity + Temperature + NO3 + Age
##
## Df AIC F Pr(>F)
## - Age 1 230.67 1.8337 0.060 .
## - NO3 1 231.25 2.3756 0.015 *
## - Salinity 1 231.62 2.7195 0.015 *
## - TOC 1 231.69 2.7842 0.015 *
## - PO4 1 232.29 3.3534 0.005 **
## - Temperature 1 232.55 3.5947 0.005 **
## - NO2 1 234.17 5.1462 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + NH4 1 231.06 1.4530 0.105
## + O2 1 231.11 1.4141 0.125
## + SPM 1 231.12 1.4041 0.185
## + HSI 1 231.44 1.1110 0.230
## + pH 1 231.38 1.1634 0.290
## + FillLevel 1 231.93 0.6668 0.840
## + SSI 1 231.99 0.6167 0.865
## + FCF 1 232.04 0.5713 0.920
## + GSI 1 232.08 0.5327 0.980
##
## [1] "Ordisep GC"
## Df AIC F Pr(>F)
## + NO2 1 240.77 6.1761 0.005 **
## + PO4 1 236.46 6.3245 0.005 **
## + TOC 1 235.16 3.2023 0.005 **
## + Salinity 1 233.51 3.5185 0.005 **
## + Temperature 1 231.72 3.6029 0.005 **
## + NO3 1 230.67 2.8494 0.010 **
## + Age 1 230.67 1.8337 0.035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Removing O2 for Collinearity"
## [1] "VIF"
## NO2 Age TOC Salinity PO4 NO3
## 1.454990 1.725781 1.774548 1.851783 2.344520 4.051192
## Temperature
## 4.052744## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## [1] "Salinity" "Temperature" "NH4" "NO2" "NO3"
## [6] "O2" "PO4" "pH" "SPM" "TOC"
## [11] "fSeason" "fLOC"
## [1] "No Waterfiler Samples removed ;ps_WF"
## [1] "NAs in TraitData"
##
## FALSE
## 190
## [1] Temperature Salinity NH4 NO2 NO3 O2
## [7] PO4 pH TOC SPM
## <0 rows> (or 0-length row.names)
## [1] "Smallest Group for avgdist"
## [1] 3576
## [1] "PERMANOVA WF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = avgdist_abund_table ~ ., data = df_Include)
## Df SumOfSqs R2 F Pr(>F)
## Temperature 1 0.8889 0.20631 8.5455 0.001 ***
## Salinity 1 0.8819 0.20467 8.4778 0.001 ***
## NH4 1 0.2410 0.05593 2.3166 0.028 *
## NO2 1 0.1149 0.02667 1.1046 0.366
## NO3 1 0.2782 0.06456 2.6741 0.011 *
## O2 1 0.2057 0.04774 1.9775 0.042 *
## PO4 1 0.2274 0.05279 2.1865 0.029 *
## pH 1 0.3006 0.06977 2.8898 0.013 *
## TOC 1 0.1443 0.03348 1.3869 0.187
## SPM 1 0.1937 0.04495 1.8620 0.085 .
## Residual 8 0.8322 0.19314
## Total 18 4.3087 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "VIF"
## Salinity PO4 pH TOC NO3 SPM
## 7.176865 9.171618 9.915326 24.715622 47.988849 53.281803
## NO2 NH4 Temperature O2
## 58.507287 108.432044 254.059907 278.295693
##
## Start: avgdist_abund_table ~ 1
##
## Df AIC F Pr(>F)
## + pH 1 25.785 5.0481 0.005 **
## + Salinity 1 26.306 4.4512 0.005 **
## + Temperature 1 26.335 4.4188 0.005 **
## + NO3 1 26.665 4.0497 0.005 **
## + O2 1 27.159 3.5098 0.005 **
## + TOC 1 28.332 2.2816 0.025 *
## + SPM 1 28.572 2.0399 0.045 *
## + PO4 1 29.051 1.5652 0.095 .
## + NO2 1 29.343 1.2821 0.175
## + NH4 1 29.669 0.9711 0.440
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ pH
##
## Df AIC F Pr(>F)
## - pH 1 28.725 5.0481 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Temperature 1 21.121 6.7209 0.005 **
## + O2 1 21.776 5.9512 0.005 **
## + NO3 1 22.398 5.2446 0.005 **
## + TOC 1 24.965 2.5599 0.015 *
## + Salinity 1 25.726 1.8311 0.055 .
## + PO4 1 25.908 1.6608 0.090 .
## + SPM 1 26.143 1.4436 0.120
## + NH4 1 26.759 0.8874 0.510
## + NO2 1 27.094 0.5923 0.940
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ pH + Temperature
##
## Df AIC F Pr(>F)
## - Temperature 1 25.785 6.7209 0.005 **
## - pH 1 26.335 7.3885 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + PO4 1 20.503 2.2162 0.005 **
## + NO3 1 20.669 2.0664 0.010 **
## + Salinity 1 21.280 1.5260 0.055 .
## + TOC 1 21.752 1.1207 0.270
## + NH4 1 21.933 0.9685 0.560
## + SPM 1 22.003 0.9096 0.565
## + NO2 1 22.135 0.7989 0.755
## + O2 1 22.275 0.6835 0.905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ pH + Temperature + PO4
##
## Df AIC F Pr(>F)
## - PO4 1 21.121 2.2162 0.030 *
## - Temperature 1 25.908 7.1490 0.005 **
## - pH 1 26.208 7.5014 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Salinity 1 20.131 1.8615 0.015 *
## + TOC 1 20.584 1.4878 0.070 .
## + NO3 1 20.690 1.4019 0.145
## + SPM 1 21.001 1.1515 0.310
## + NO2 1 21.152 1.0314 0.400
## + NH4 1 21.364 0.8649 0.620
## + O2 1 21.580 0.6970 0.860
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: avgdist_abund_table ~ pH + Temperature + PO4 + Salinity
##
## Df AIC F Pr(>F)
## - Salinity 1 20.503 1.8615 0.065 .
## - PO4 1 21.280 2.5240 0.010 **
## - pH 1 22.071 3.2262 0.005 **
## - Temperature 1 25.327 6.4457 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + NO3 1 20.052 1.5030 0.060 .
## + NO2 1 20.569 1.1143 0.300
## + SPM 1 20.556 1.1233 0.305
## + TOC 1 20.637 1.0637 0.380
## + NH4 1 20.753 0.9780 0.425
## + O2 1 21.073 0.7443 0.775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] "Ordisep WF"
## Df AIC F Pr(>F)
## + pH 1 25.785 5.0481 0.005 **
## + Temperature 1 21.121 6.7209 0.005 **
## + PO4 1 20.503 2.2162 0.005 **
## + Salinity 1 20.131 1.8615 0.015 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Removing O2 for Collinearity"
## [1] "VIF"
## Temperature PO4 pH Salinity
## 1.197473 1.199426 3.014100 3.091885## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
pslist$RDA_list <- RDA_list
print(sort(vif.cca(RDA_list$dbRDA_results_OE)))## FillLevel SPM HSI PO4 TOC GSI
## 1.135031 1.471549 1.778025 2.008939 2.489322 2.839159
## pH Salinity Temperature NO2 NO3 NH4
## 3.096458 5.051252 6.326722 6.891286 6.903946 7.898809print(sort(vif.cca(RDA_list$dbRDA_results_GC)))## NO2 Age TOC Salinity PO4 NO3
## 1.454990 1.725781 1.774548 1.851783 2.344520 4.051192
## Temperature
## 4.052744#Export EnvFit Results
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
Datasetname <- sub("ps_", "", Dataset)
EnvFit_results<- as.data.frame(
cbind(round(pslist$RDA_list[[paste("EnvFit_results", Datasetname, sep="_")]]$vectors$r, digits=3),
round(pslist$RDA_list[[paste("EnvFit_results", Datasetname, sep="_")]]$vectors$pvals, digits=3))
) %>%
setNames(c("R2", "P"))
write.csv2(EnvFit_results, paste0(file.path(path_Output_16S, paste(paste(save_name,
paste("EnvFit_results", sep=""), Datasetname, sep="_"),".csv", sep=""))))
}6.3.5 Mantel Test
https://jkzorz.github.io/2019/07/08/mantel-test.html Mantel tests are correlation tests that determine the correlation between two matrices (rather than two variables). When using the test for microbial ecology, the matrices are often distance/dissimilarity matrices with corresponding positions (i.e. samples in the same order in both matrices). In order to calculate the correlation, the matrix values of both matrices are ‘unfolded’ into long column vectors, which are then used to determine correlation. Permutations of one matrix are used to determine significance. Learn more about Mantel tests here.
#########################
#Fish_OTU vs Environment#
#########################
#based in the avgdist matrix
Mantel_list <- list()
for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis ))]) {
Mantel_list_length <- length(Mantel_list)
Datasetname <- sub("_Avg_dist_matrix", "", Dataset)
ps <- pslist[[paste("ps", Datasetname, sep="_")]]
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
metadata <- as.data.frame(sample_data(ps_Filter))
metadata <- metadata[,names(metadata) %in% c("SampleID", "Temperature", "Salinity", "O2", "SPM", "PO4", "NH4", "NO3", "NO2")]
dist_tibble <- pslist$BrayCurtis[[Dataset]] %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample)
mantel_function = function(species, parameter){
rare_dist_matrix <- dist_tibble %>%
left_join(., metadata, by = c("sample" = "SampleID"))%>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, name, value) %>%
pivot_wider(names_from = name, values_from = value) %>%
select(-sample) %>%
as.dist()
env_dist <- dist_tibble %>%
left_join(., metadata, by = c("sample" = "SampleID")) %>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, paste0(parameter,".x")) %>%
unique() %>%
column_to_rownames("sample") %>%
dist(method="euclidean")
abund_temp = mantel(rare_dist_matrix, env_dist, method = "spearman", permutations = 9999, na.rm = TRUE)
cbind(species=species, env=parameter, mantel=abund_temp$statistic, sig=abund_temp$signif)
}
Env.out = data.frame(species=NA, env=NA, mantel=NA, sig=NA)
species=Datasetname
envs=c("Temperature", "Salinity", "O2", "SPM", "PO4", "NH4", "NO3", "NO2")
for(Dat in species){
for(params in envs)
Env.out <- rbind(Env.out, mantel_function(species = Dat, parameter = params))
}
Mantel_list[[Mantel_list_length+1]] <- Env.out
names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "Env", sep="_")
}
########################
#Fish_OTU vs Physiology#
########################
#based in the avgdist matrix
for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis))]) {
Mantel_list_length <- length(Mantel_list)
Datasetname <- sub("_Avg_dist_matrix", "", Dataset)
ps <- pslist[[paste("ps", Datasetname, sep="_")]]
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
Phys <- c(
"Length", "Age", "FCF", "GSI", "HSI", "SSI", "FillLevel")
'%!in%' <- function(x,y)!('%in%'(x,y))
df_Traits <- data.frame(sample_data(ps_Filter))
df_Traits <- df_Traits[names(df_Traits) %in% Phys]
print("NAs in TraitData")
print(table(is.na(df_Traits)))
print((df_Traits[!complete.cases(df_Traits), ]))
df_Traits <- df_Traits[complete.cases(df_Traits), ]
df_Traits <- decostand(df_Traits, method='standardize')
ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
rownames(df_Traits[complete.cases(df_Traits), ])], ps_Filter)
##################################################
#Determine which variables are useful for the plot#
##################################################
abund_table <- as.data.frame(otu_table(ps_Filter, taxa_are_rows = FALSE))
smallest_group <- min(rowSums(otu_table(ps_Filter)))
avgdist_abund_table <- avgdist(abund_table, dmethod = "bray", sample = smallest_group)
dist_tibble <- avgdist_abund_table %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample)
metadata <- as.data.frame(sample_data(ps_Filter))
metadata <- metadata[,names(metadata) %in% c("SampleID", Phys)]
mantel_function = function(species, parameter){
rare_dist_matrix <- dist_tibble %>%
left_join(., metadata, by = c("sample" = "SampleID"))%>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, name, value) %>%
pivot_wider(names_from = name, values_from = value) %>%
select(-sample) %>%
as.dist()
env_dist <- dist_tibble %>%
left_join(., metadata, by = c("sample" = "SampleID")) %>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, paste0(parameter,".x")) %>%
unique() %>%
column_to_rownames("sample") %>%
dist(method="euclidean")
abund_temp = mantel(rare_dist_matrix, env_dist, method = "spearman", permutations = 9999, na.rm = TRUE)
cbind(species=species, env=parameter, mantel=abund_temp$statistic, sig=abund_temp$signif)
}
Phys.out = data.frame(species=NA, env=NA, mantel=NA, sig=NA)
species=Datasetname
for(Dat in species){
for(params in Phys)
Phys.out <- rbind(Phys.out, mantel_function(species = Dat, parameter = params))
}
Mantel_list[[Mantel_list_length+1]] <- Phys.out
names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "Phys", sep="_")
} ## [1] "NAs in TraitData"
##
## FALSE
## 966
## [1] Length FillLevel Age HSI SSI GSI FCF
## <0 rows> (or 0-length row.names)
## [1] "NAs in TraitData"
##
## FALSE
## 602
## [1] Length FillLevel Age HSI SSI GSI FCF
## <0 rows> (or 0-length row.names)#########################
#Water_OTU vs Environment#
#########################
#based in the avgdist matrix
for (Dataset in names(pslist)[grepl("ps_WF", names(pslist))]) {
Mantel_list_length <- length(Mantel_list)
Datasetname <- sub("ps_", "", Dataset)
ps <- pslist[[paste("ps", Datasetname, sep="_")]]
ps_Filter <- ps
print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
print("Waterfilters will not be filtered")
Envs <- c(
"Temperature", "Salinity", "O2", "SPM", "PO4", "NH4", "NO3", "NO2")
'%!in%' <- function(x,y)!('%in%'(x,y))
df_Traits <- data.frame(sample_data(ps_Filter))
df_Traits <- df_Traits[names(df_Traits) %in% Envs]
print("NAs in TraitData")
print(table(is.na(df_Traits)))
print((df_Traits[!complete.cases(df_Traits), ]))
df_Traits <- df_Traits[complete.cases(df_Traits), ]
df_Traits <- decostand(df_Traits, method='standardize')
ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
rownames(df_Traits[complete.cases(df_Traits), ])], ps_Filter)
##################################################
#Determine which variables are useful for the plot#
##################################################
abund_table <- as.data.frame(otu_table(ps_Filter, taxa_are_rows = FALSE))
smallest_group <- min(rowSums(otu_table(ps_Filter)))
print("Smallest group for avgdist")
print(smallest_group)
avgdist_abund_table <- avgdist(abund_table, dmethod = "bray", sample = smallest_group)
dist_tibble <- avgdist_abund_table %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample)
metadata <- as.data.frame(sample_data(ps_Filter))
metadata <- metadata[,names(metadata) %in% c("SampleID", Envs)]
mantel_function = function(species, parameter){
rare_dist_matrix <- dist_tibble %>%
left_join(., metadata, by = c("sample" = "SampleID"))%>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, name, value) %>%
pivot_wider(names_from = name, values_from = value) %>%
select(-sample) %>%
as.dist()
env_dist <- dist_tibble %>%
left_join(., metadata, by = c("sample" = "SampleID")) %>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, paste0(parameter,".x")) %>%
unique() %>%
column_to_rownames("sample") %>%
dist(method="euclidean")
abund_temp = mantel(rare_dist_matrix, env_dist, method = "spearman", permutations = 9999, na.rm = TRUE)
cbind(species=species, env=parameter, mantel=abund_temp$statistic, sig=abund_temp$signif)
}
Env.out = data.frame(species=NA, env=NA, mantel=NA, sig=NA)
species=Datasetname
envs=Envs
for(Dat in species){
for(params in envs)
Env.out <- rbind(Env.out, mantel_function(species = Dat, parameter = params))
}
Mantel_list[[Mantel_list_length+1]] <- Env.out
names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "Env", sep="_")
} ## [1] "Samples removed for low frequency below 5000 seqs in;ps_WF"
## WFSU21SSFL WFWI22MGEB
## 3 6
## [1] "Waterfilters will not be filtered"
## [1] "NAs in TraitData"
##
## FALSE
## 152
## [1] Temperature Salinity NH4 NO2 NO3 O2
## [7] PO4 SPM
## <0 rows> (or 0-length row.names)
## [1] "Smallest group for avgdist"
## [1] 3576#######################
#Fish_ASV vs Water_ASV#
#######################
#-> Problem No Waterfilters for Autumn and SPTW
#Distance Matrix Fish
#Distance Matrix Water
#->Elongate Water Matrix to fit the fish Matrix
# for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis ))]) {
#
# Mantel_list_length <- length(Mantel_list)
# Datasetname <- sub("_Avg_dist_matrix", "", Dataset)
#
# ps_Fish <- pslist[[paste("ps", Datasetname, sep="_")]]
# ps_Filter_Fish <- subset_samples(ps_Fish, sample_sums(ps_Fish) > 5000)
#
# ps_WF <- pslist[[paste("ps", "WF", sep="_")]]
# abund_table_WF <- data.frame(otu_table(ps_WF))
#
# ###########################
# #Avg_dist for Fish samples#
# ###########################
# abund_table_Fish <- as.data.frame(otu_table(ps_Filter_Fish, taxa_are_rows = FALSE))
# print(Datasetname)
# print("Smalles Group AVG dist")
# smallest_group_Fish <- min(rowSums(otu_table(ps_Filter_Fish)))
# print(smallest_group_Fish)
# avgdist_abund_table_Fish <- avgdist(abund_table_Fish, dmethod = "bray", sample = smallest_group_Fish)
#
# dist_tibble_Fish <- avgdist_abund_table_Fish %>%
# as.matrix() %>%
# as_tibble(rownames="sample") %>%
# pivot_longer(-sample)
#
# ################################################
# #Avg_fist for artificially matched Watersamples#
# ################################################
# Names_Fish <- substring(sample_names(ps_Filter_Fish), first = 3)
# WF_mantel<- list()
# for (i in seq_along(Names_Fish)) {
# NAME <- Names_Fish[i]
# NAME2 <- gsub("\\d+$", "", NAME)
# NAME3 <- substr(NAME2 , 1, nchar(NAME2 ) - 2)
# NAME4 <- substr(NAME3 , 1, nchar(NAME3 ) - 2)
#
# if (paste("WF", NAME3, "EB", sep="") %in% rownames(abund_table_WF)) {
# WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "EB",
# sep=""),]
#
# } else if (paste("WF", NAME3, "FL", sep="") %in% rownames(abund_table_WF)){
# WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "FL",
# sep=""),]
#
# } else if (paste("WF", NAME4, "TWEB", sep="") %in% rownames(abund_table_WF)){
# WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "TWEB",
# sep=""),]
#
# } else if (paste("WF", NAME4, "MLFL", sep="") %in% rownames(abund_table_WF)){
# WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "MLFL",
# sep=""),]
#
# } else {print(NA)}
#
# names(WF_mantel)[[i]] <- paste("WF", NAME, sep="")
# }
#
# WF_mantel_df <- do.call(rbind, WF_mantel)
#
# print("Smallest Group AVG dist Waterfilter")
# smallest_group_WF <- min(rowSums(otu_table(ps_WF)))
# print(smallest_group_WF)
# avgdist_abund_table_WF <- avgdist(WF_mantel_df , dmethod = "bray", sample = smallest_group_WF)
#
# #####################
# #Run the Mantel Test#
# #####################
# dist_tibble_Fish <- avgdist_abund_table_Fish %>%
# as.matrix() %>%
# as_tibble(rownames="sample") %>%
# pivot_longer(-sample)
#
# dist_tibble_WF <- avgdist_abund_table_WF %>%
# as.matrix() %>%
# as_tibble(rownames="sample") %>%
# pivot_longer(-sample)
#
# rare_dist_matrix_Fish <- dist_tibble_Fish %>%
# left_join(., metadata, by = c("sample" = "SampleID"))%>%
# left_join(., metadata, by = c("name" = "SampleID")) %>%
# select(sample, name, value) %>%
# pivot_wider(names_from = name, values_from = value) %>%
# select(-sample) %>%
# as.dist()
#
# rare_dist_matrix_WF <- dist_tibble_WF %>%
# left_join(., metadata, by = c("sample" = "SampleID"))%>%
# left_join(., metadata, by = c("name" = "SampleID")) %>%
# select(sample, name, value) %>%
# pivot_wider(names_from = name, values_from = value) %>%
# select(-sample) %>%
# as.dist()
#
# abund_temp <-
# mantel(rare_dist_matrix_Fish, rare_dist_matrix_WF, method = "spearman", permutations = 9999, na.rm = TRUE)
# Env.out = data.frame(species=Datasetname, env="WF", mantel=abund_temp$statistic, sig=abund_temp$signif)
#
# Mantel_list[[Mantel_list_length+1]] <- Env.out
# names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "WF", sep="_")
# }
##########################################################################
#Without Autumn Samples as there is a large lack in matching watersamples#
##########################################################################
for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis ))]) {
Mantel_list_length <- length(Mantel_list)
Datasetname <- sub("_Avg_dist_matrix", "", Dataset)
ps_Fish <- pslist[[paste("ps", Datasetname, sep="_")]]
ps_Filter_Fish <- subset_samples(ps_Fish, sample_sums(ps_Fish) > 5000)
ps_Filter_Fish <- prune_samples(!grepl("AU21BB|AU21MG|AU21ML|AU21SS|SU21TW|WI22TW|SP22TW", sample_names(ps_Filter_Fish)), ps_Filter_Fish)
ps_WF <- pslist[[paste("ps", "WF", sep="_")]]
abund_table_WF <- data.frame(otu_table(ps_WF))
###########################
#Avg_dist for Fish samples#
###########################
abund_table_Fish <- as.data.frame(otu_table(ps_Filter_Fish, taxa_are_rows = FALSE))
print(Datasetname)
print("Smalles Group AVG dist")
smallest_group_Fish <- min(rowSums(otu_table(ps_Filter_Fish)))
print(smallest_group_Fish)
avgdist_abund_table_Fish <- avgdist(abund_table_Fish, dmethod = "bray", sample = smallest_group_Fish)
dist_tibble_Fish <- avgdist_abund_table_Fish %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample)
################################################
#Avg_fist for artificially matched Watersamples#
################################################
#we lack some watersamples matching with all fish data, as the bacterioplankton appears rather homogenous compared to the fish, we artificially match water samples from the closest time and space with the fish data
Names_Fish <- substring(sample_names(ps_Filter_Fish), first = 3)
WF_mantel<- list()
for (i in seq_along(Names_Fish)) {
NAME <- Names_Fish[i]
NAME2 <- gsub("\\d+$", "", NAME)
NAME3 <- substr(NAME2 , 1, nchar(NAME2 ) - 2)
NAME4 <- substr(NAME3 , 1, nchar(NAME3 ) - 2)
if (paste("WF", NAME3, "EB", sep="") %in% rownames(abund_table_WF)) {
WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "EB",
sep=""),]
} else if (paste("WF", NAME3, "FL", sep="") %in% rownames(abund_table_WF)){
WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "FL",
sep=""),]
} else if (paste("WF", NAME4, "TWEB", sep="") %in% rownames(abund_table_WF)){
WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "TWEB",
sep=""),]
} else if (paste("WF", NAME4, "MLFL", sep="") %in% rownames(abund_table_WF)){
WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "MLFL",
sep=""),]
} else {print(NA)}
names(WF_mantel)[[i]] <- paste("WF", NAME, sep="")
}
WF_mantel_df <- do.call(rbind, WF_mantel)
print("Smallest Group AVG dist Waterfilter")
smallest_group_WF <- min(rowSums(otu_table(ps_WF)))
print(smallest_group_WF)
avgdist_abund_table_WF <- avgdist(WF_mantel_df , dmethod = "bray", sample = smallest_group_WF)
#####################
#Run the Mantel Test#
#####################
dist_tibble_Fish <- avgdist_abund_table_Fish %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample)
dist_tibble_WF <- avgdist_abund_table_WF %>%
as.matrix() %>%
as_tibble(rownames="sample") %>%
pivot_longer(-sample)
rare_dist_matrix_Fish <- dist_tibble_Fish %>%
left_join(., metadata, by = c("sample" = "SampleID"))%>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, name, value) %>%
pivot_wider(names_from = name, values_from = value) %>%
select(-sample) %>%
as.dist()
rare_dist_matrix_WF <- dist_tibble_WF %>%
left_join(., metadata, by = c("sample" = "SampleID"))%>%
left_join(., metadata, by = c("name" = "SampleID")) %>%
select(sample, name, value) %>%
pivot_wider(names_from = name, values_from = value) %>%
select(-sample) %>%
as.dist()
abund_temp <-
mantel(rare_dist_matrix_Fish, rare_dist_matrix_WF, method = "spearman", permutations = 9999, na.rm = TRUE)
Env.out = data.frame(species=Datasetname, env="WF", mantel=abund_temp$statistic, sig=abund_temp$signif)
Mantel_list[[Mantel_list_length+1]] <- Env.out
names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "WF", "NoAutumn",sep="_")
}## [1] "OE"
## [1] "Smalles Group AVG dist"
## [1] 9770
## [1] "Smallest Group AVG dist Waterfilter"
## [1] 3576
## [1] "GC"
## [1] "Smalles Group AVG dist"
## [1] 6441
## [1] "Smallest Group AVG dist Waterfilter"
## [1] 3576pslist[["Mantel_list"]] <- Mantel_list
# > Mantel_list$`OE_WF_MG-BB`
# species env mantel sig
# 1 OE WF 0.5307317 1e-04
# > Mantel_list$`GC_WF_MG-BB`
# species env mantel sig
# 1 GC WF 0.3847558 0.0149
# Mantel_list$`OE_WF_SS-ML`
# species env mantel sig
# 1 OE WF 0.2997895 1e-04
# Mantel_list$`GC_WF_SS-ML`
# species env mantel sig
# 1 GC WF 0.07524223 0.1884
#####################
#Plot Mantel Results#
#####################
Mantel_list_clean <- lapply(Mantel_list, function(x) {
x[complete.cases(x), ]
})
Mantel_out <- as.data.frame(do.call(rbind, Mantel_list_clean))
rownames(Mantel_out) <- NULL
Mantel_out$sig <- as.numeric(Mantel_out$sig)
Mantel_out$mantel <- as.numeric(Mantel_out$mantel)
Mantel_out_sig <- Mantel_out %>%
filter(sig < 0.05)
vars <- c("Salinity", "Temperature", "O2","SPM", "NO3", "NH4", "NO2", "PO4",
"FCF", "GSI", "HSI", "Length", "Age", "WF")
Envs <- c("Salinity", "Temperature", "O2","SPM", "NO3", "NH4", "NO2", "PO4")
Phy <- c("Length", "Age", "FCF","HSI", "GSI")
Mantel_out_sig <- Mantel_out_sig %>%
mutate(Species= factor(species, levels = c("OE", "GC", "WF"))) %>%
mutate(Env=factor(env, levels = vars)) %>%
mutate(kind = case_when(
env %in% Envs ~ "Environment",
env %in% c("WF") ~ "ASV",
TRUE ~ "Physiology"))
Mantel_plot <- Mantel_out_sig %>%
ggplot(., aes(x=Env, y=Species, fill=mantel)) +
geom_tile() +
scale_fill_gradient2(
low = "#FFFFFF",
high = "#006000",
limit = c(0,1)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, size = 15, face ="bold"),
#axis.title.x.bottom = element_text(size=15,face = "bold"),
#axis.title.y.left = element_text(size=15,face = "bold"),
strip.text.x = element_text(size = 15, face = "bold"),
strip.text.x.bottom = element_text(size = 15,face = "bold"),
strip.text.y.left = element_text(size = 15,face = "bold"),
axis.title.x.bottom = element_blank(),
axis.title.y.left = element_blank(),
axis.text.x.bottom = element_text(face = "bold", angle = 45, hjust = 1),
axis.text.y.left = element_text(face = "bold", size = 15),
legend.title = element_text(size = 15, face ="bold"),
legend.text = element_text(size = 15,face = "bold"),
plot.title = element_text(size = 15, face = "bold"),
plot.subtitle = element_text(size = 15),
plot.caption = element_text(size = 15)) +
labs(#title = "Correlation relationship between avg.dist Bray-Curtis matrices and varibles", y = "Sample kind",
fill="corr") +
labs(fill = "Mantel's r") +
#theme(legend.box.background = element_rect(color = "black"))
theme(legend.position="right", legend.key.width = unit(1, "cm"))
Mantel_plot <- Mantel_plot + geom_text(vjust=0.76,hjust=0.49, aes(label = paste0(c("","*")[(abs(sig) <= .05)+1],
c("","*")[(abs(sig) <= .001)+1], c("","*")[(abs(sig) <= .001)+1])))
kindOrder <- c("Environment", "Physiology", "ASV")
Mantel_plot <- Mantel_plot + facet_grid(. ~ factor(kind, levels = kindOrder), drop = TRUE, scale = "free", space = "free_x")
Mantel_plot <- cowplot::plot_grid(Mantel_plot, labels = c(""), ncol = 1)
Mantel_plot <- plot_grid(Mantel_plot, ncol = 1, rel_heights = c(0.02, 0.05, 0.98))
plot(Mantel_plot)
ggsave(Mantel_plot, filename = paste(paste("Mantel_test",
sep="_"),".png", sep=""), path = pathPlots , device='png', dpi=300, width = 7,height = 3)6.4 Export ASV & Core Dataframe
#############################################
#Create Bacterial Counts and Reseq Dataframe#
#############################################
SSU_Data_list <- list()
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
require(plyr)
require(dplyr)
require(phyloseq)
Datasetname <- sub("ps_", "", Dataset)
ps <- pslist[[Dataset]]
SSU_Data_list_length <- length(SSU_Data_list)
TAX <- as.data.frame(tax_table(ps))
OTU <- as.data.frame(t(otu_table(ps)))
REFSEQ <- refseq(ps)
REFSEQ <- stack(as.character(REFSEQ, use.names=TRUE))
colnames(REFSEQ)[colnames(REFSEQ) == "ind"] <- "ASV"
#################################
#Creating AverageSums per Season#
#################################
SeasonSums <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
as.data.frame() %>%
rownames_to_column(var = "SampleID") %>%
dplyr::left_join(SAMDF16S, by = c("SampleID" = "SampleID")) %>%
dplyr::group_by(Season) %>%
dplyr::summarise(dplyr::across(rownames(tax_table(ps)), mean, na.rm = TRUE)) %>%
t() %>%
as.data.frame() %>%
`colnames<-`(.[1, ]) %>%
.[-1, ] %>%
setNames(paste0("avg_", colnames(.))) %>%
mutate_all(as.numeric) %>%
rownames_to_column(var="ASV")
SSU_Data <- TAX %>%
rownames_to_column(var = "ASV") %>%
left_join( #Add relative ASVmeans
data.frame(t(phyloseq::otu_table(ps %>%
transform_sample_counts(function(x) {x/sum(x)*100})))) %>%
dplyr::mutate(ASVmeans = rowMeans(.)) %>%
dplyr::mutate(ASV = rownames(.)) %>%
dplyr::select(ASV, ASVmeans)
) %>%
left_join( #Add Sums by Season
SeasonSums
) %>%
left_join( #Add relative ASV data per sample
data.frame(t(phyloseq::otu_table(ps %>%
transform_sample_counts(function(x) {x/sum(x)*100})))) %>%
rownames_to_column(var = "ASV")
) %>%
left_join(REFSEQ
) %>%
dplyr::arrange(desc(ASVmeans))
rownames(SSU_Data) <- SSU_Data$ASV
write.csv2(SSU_Data, paste0(file.path(path_Output_16S, paste(paste(save_name,
paste("SSU_ALL", sep=""), Date, Datasetname, sep="_"),".csv", sep=""))))
###################
#Extract Core-Taxa#
###################
core_taxa_standard <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
microbiome::core_members(., detection = 0.0001, prevalence = 90/100)
# ps_rel <- microbiome::transform(pslist[[paste("ps_", Species, sep="")]], "compositional")
# #Initial core computing
# #Setting a detection threshold of 0.0001% (within samples) and prevelance threshold of
#80% (across samples).
# core_taxa_standard <- microbiome::core_members(ps_rel, detection = 0.0001, prevalence = 90/100)
write.csv2(SSU_Data[SSU_Data$ASV%in% core_taxa_standard,],
file.path(path_Output_16S,paste(paste(save_name,"16S_merge_Filter_ASV_besttax_Core",
Datasetname, Date, sep="_"), ".csv", sep="")))
SSU_Data_list[[SSU_Data_list_length+1]] <- SSU_Data
names(SSU_Data_list)[[SSU_Data_list_length+1]] <- Datasetname
}#-
Save Data
###########
#Save Data#
###########
saveRDS(pslist, file.path(path_Output_16S, paste(paste(save_name, "ps_16S_merge_Filter_ASV_besttax", Date, sep="_"), ".rds", sep="")))